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Solution  of  the  Multidimensional  Compressible 
Havier- Stokes  Equations  by  a Generalized. 
Implicit  Method 


SUMMARY 


In  an  effort  to  exploit  the  favorable  stability  properties  of  implicit  methods 
and  thereby  increase  computational  efficiency  by  taking  large  time  steps,  an  implicit 
finite-difference  method  for  the  multidimensional  Navier-Stokes  equations  is  pre- 
sented. The  method  consists  of  a generalized  implicit  scheme  which  has  been  linear- 
ized by  Taylor  expansion  about  the  solution  at  the  known  time  level  to  produce  a 
set  of  coupled  linear  difference  equations  which  are  valid  for  a given  time  step. 

To  solve  these  difference  equations,  the  Douglas-Gunn  procedure  for  generating 
alternating-direction  implicit  (ADI)  schemes  as  perturbations  of  fundamental  impli- 
cit difference  schemes  is  employed.  The  resulting  sequence  of  narrow  block-banded 
systems  can  be  solved  efficiently  by  standard  block-elimination  methods.  The  method 
is  a one-step  method,  as  opposed  to  a predictor-corrector  method,  and  requires  no 
iteration  to  compute  the  solution  for  a single  time  step.  The  use  of  both  second 
and  fourth  order  spatial  differencing  is  discussed.  Test  calculations  are  presented 
for  a three-dimensional  application  to  subsonic  flow  in  a straight  duct  with  rec- 
tangular cross  section.  Stability  is  demonstrated  for  time  steps  which  are  orders 
of  magnitude  larger  than  the  maximum  allowable  time  step  for  conditionally  stable 
methods  as  determined  by  the  well  known  CFL  condition.  The  computational  effort  per 
time  step  is  discussed  and  is  very  approximately  only  twice  that  of  most  explicit 
methods.  The  accuracy  of  computed  solutions  is  examined  by  mesh  refinement  and 
comparison  with  other  analytical  and  experimental  results.  Finally,  some  test  cal- 
culations for  turbulent  flow  were  made  using  a simple  turbulence  model  consisting 
of  an  eddy  viscosity  and  specified  mixing  length.  The  results  of  these  calculations 
are  discussed. 
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INTRODUCTION 


One  of  the  major  obstacles  to  the  routine  numerical  solution  of  the 
multidimensional  compressible  Navier-Stokes  equations  is  the  large  amount  of  computer 
time  generally  required,  and  consequently,  efficient  computational  methods  are  highly 
desirable  in  this  instance.  Most  previous  methods  for  solving  the  compressible 
Navier-Stokes  equations  have  been  based  on  explicit  difference  schemes  for  the 
unsteady  form  of  the  governing  equations  and  are  subject  to  one  or  more  stability 
restrictions  on  the  size  of  the  time  step  relative  to  the  spatial  mesh  size.  These 
stability  limits  usually  correspond  to  the  well  known  Courant-Friedrichs-Lewy  (CFL) 
condition  and,  in  some  schemes,  to  an  additional  stability  condition  arising  from 
viscous  terms.  In  one  dimension,  the  CFL  condition  is  At  < Ax/(|u|  + c),  and  the 
viscous  stability  condition  is  At  < (ax)2/2v,  where  At  is  the  time  step.  Ax  is  the 
mesh  size,  u is  velocity,  c is  the  speed  of  sound,  and  v is  kinematic  viscosity. 

These  stability  restrictions  can  lower  computational  efficiency  by  imposing  a smaller 
time  step  than  would  otherwise  be  desirable.  Thus,  a key  disadvantage  of  conditionally 
stable  methods  is  that  the  maximum  time  step  is  fixed  by  the  spatial  mesh  size 
rather  than  the  physical  time  dependence  or  the  desired  temporal  accuracy.  In  con- 
trast to  most  explicit  methods,  implicit  methods  tend  to  be  stable  for  large  time 
steps  and  hence  offer  the  prospect  of  substantial  increases  in  computational 
efficiency,  provided  of  course  that  large  time  steps  are  acceptable  for  the  physical 
problem  of  interest  and  that  the  computational  effort  per  time  step  is  competitive 
with  that  of  explicit  methods.  In  an  effort  to  exploit  these  potentially  favorable 
stability  properties,  an  efficient  implicit  method  based  on  alternating-direction 
differencing  techniques  was  developed  and  is  presented  herein. 


Effect  of  Stability  Conditions 

The  potential  of  implicit  methods  for  increased  overall  efficiency  arises 
whenever  the  CFL,  viscous,  orofcher  stability  conditions  become  restrictive  in  the 
sense  that  the  rate  of  change  of  physical  processes  permits  a larger  step  than  does 
the  numerical  scheme.  Since  the  severity  of  these  stability  conditions  varies  from 
problem  to  problem  and  with  the  choice  of  grid  size,  no  single  factor  can  be  quoted 
for  overall  relative  efficiency.  Nevertheless,  certain  guidelines  can  be  established 
from  a consideration  of  the  effect  various  factors  have  on  the  stability  conditions. 
One  of  the  adverse  consequences  of  stability  restrictions  becomes  apparent  when  a 
coarse-mesh  solution  is  recomputed  with  a finer  spatial  mesh  to  obtain  greater 
spatial  accuracy.  Assumirg  the  maximum  allowable  time-step  is  taken  with  a condi- 
tionally stable  scheme,  the  time-step  must  be  smaller  with  the  finer  mesh  even 
though  the  physical  time  dependence  is  exactly  the  same.  In  these  same  circumstances, 
implicit  methods  not  subject  to  stability  restrictions  can  take  the  same  time  step 
with  both  spatial  meshes,  and  consequently,  the  computational  effort  increases  only 
linearly  with  the  number  of  spatial  grid  points,  as  opposed  to  a quadratic  increase 


2 


R75-9U363-15 


for  a CFL-limited  calculation  and  a cubic  increase  for  a viscous-limited  calculation. 
For  lower  Mach  number  flows,  the  CFL  condition  eventually  becomes  restrictive  regard- 
less of  mesh  spacing,  and  thus  stable  implicit  methods  are  well  suited  for  problems 
in  the  low  Mach  number  regime.  For  all  Mach  numbers,  both  the  CFL  and  viscous 
stability  conditions  eventually  become  restrictive  for  sufficiently  small  mesh 
spacing.  Implicit  methods  thus  become  increasingly  attractive  when  high  spatial 
resolution  is  necessary  and  especially  when  locally  refined  meshes  are  used,  since 
the  stability  limit  is  usually  governed  by  the  smallest  mesh  spacing  in  the  field. 
This  latter  situation  is  common  when  more  than  one  length  scale  is  present,  as 
is  the  case  for  flows  which  are  largely  inviscid  but  have  thin  boundary  layers 
requiring  a locally  refined  mesh.  Similar  statements  hold  for  (time -averaged)  tur- 
bulent flows  with  viscous  sublayers . 


Steady  and  Unsteady  Applications 

Although  the  present  method  solves  the  unsteady  equations  of  motion,  the  most 
dramatic  gains  in  efficiency  are  likely  to  result  when  computing  steady  solutions 
as  the  asymptotic  limit  for  large  time  of  an  unsteady  solution.  In  this  instance, 
there  is  no  need  to  follow  the  transient  accurately,  and  large  temporal  truncation 
errors  can  therefore  be  tolerated  and  are  even  desirable  in  exchange  for  a reduction 
in  the  total  number  of  time  steps  required  to  reach  steady  conditions.  Furthermore, 
various  convergence  acceleration  techniques  which  take  advantage  of  stability  can 
be  employed.  For  unsteady  problems,  the  need  for  transient  accuracy  limits  the 
extent  to  which  the  time  step  can  be  increased  as  a result  of  improved  stability, 
and  the  relevant  question  then  becomes  how  the  stability  limit  compares  with  the 
time  step  necessary  to  follow  the  transient  accurately.  This  question  is  difficult 
to  answer  with  any  generality  since  stability  conditions  only  relate  the  time  step 
to  the  spatial  mesh  size,  and  they  do  not  directly  relate  computational  time  step 
to  the  relevant  physical  time  scales  of  the  overall  problem. 

Although  the  use  of  implicit  schemes  for  diffusive  terms  is  generally  accepted, 
it  is  less  clear  whether  implicit  schemes  for  convective  terms  have  sufficient 
transient  accuracy  for  unsteady  problems (cf.  Orszag  & Israeli,  Ref.  l) . Aside 
from  the  usual  order-of-accuracy  error  estimates,  which  are  relevant  only  for 
sufficiently  small  step  size,  indications  of  the  accuracy  of  a scheme  for  a given 
finite  step  are  valuable.  Thus,  various  authors  [e.g.,  Fromm  (Ref.  2);  Morton 
(Ref.  3)]  have  considered  the  damping  and  dispersion  characteristics  of  different 
difference  schemes  in  approximating  pure-convection  phenomena.  For  a model  problem, 
Morton  (Ref.  3)  has  pointed  out  that  while  the  Crank -Nicolson  implicit  scheme  intro- 
duces no  damping  error,  certain  explicit  schemes  perform  significantly  better  in  terms 
of  phase  error  per  time  step  as  a function  of  wavelength  than  does  the  Crank- 
Nicolson  scheme.  Furthermore,  it  can  be  shown  for  this  model  problem  that  the  phase 
error  of  the  Crank-Nicolson  scheme  becomes  progressively  worse  as  the  explicit 
stability  limit  is  exceeded,  even  though  the  scheme  itself  remains  stable.  These 
findings  suggest  that  implicit  schemes  for  convective  terms  may  not  have  sufficient 
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transient  accuracy  for  unsteady  problems  unless  they  also  satisfy  the  explicit 
stability  conditions.  This  argument  no  doubt  has  validity  in  numerous  applications 
wherein  the  physical  time  scale  of  interest  is  the  same  as  the  time  scale  governing 
stability.  However,  the  pure-convection  model  problem  does  not  explore  the  effects 
of  multiple  time  scales,  boundary  conditions,  or  local  regions  of  high  resolution, 
which  are  often  important  in  practical  applications.  Thus,  for  example,  low  Mach 
number  unsteady  flows  can  be  envisaged  wherein  the  velocity  field  is  insensitive  to 
the  detailed  behavior  of  sound  waves.  Similarly,  the  unsteadiness  may  enter 
through  periodic  boundary  conditions  having  a time  scale  longer  than  that  governing 
stability.  Further  study  in  this  area  would  be  helpful  in  assessing  the  value  of 
implicit  schemes  in  unsteady  applications. 


The  Present  Method 


Summary 

The  present  method  can  be  briefly  outlined  as  follows:  the  governing  equations 

are  replaced  by  an  implicit  time  difference  approximation,  optionally  a backward 
difference  or  Crank -Nicolson  scheme.  Terms  involving  nonlinearities  at  the  implicit 
time  level  are  linearized  by  Taylor  expansion  about  the  solution  at  the  known  time 
level,  and  spatial  difference  approximations  are  introduced.  The  result  is  a system 
of  multidimensional  coupled  (but  linear)  difference  equations  for  the  dependent 
variables  at  the  unknown  or  implicit  time  level.  To  solve  these  difference 
equations,  the  Douglas-Gunn  (Ref.  4)  procedure  for  generating  alternating-direction 
implicit  (ADI)  schemes  as  perturbations  of  fundamental  implicit  difference  schemes 
is  introduced.  This  technique  leads  to  systems  of  coupled  linear  difference  equa- 
tions having  narrow  block-banded  matrix  structures  which  can  be  solved  efficiently 
by  standard  block- elimination  methods. 

The  method  centers  around  the  use  of  a formal  linearization  technique  adapted 
for  the  integration  of  initial -value  problems.  The  linearization  technique,  which 
of  necessity  requires  an  implicit  solution  procedure,  permits  the  solution  of  coupled 
nonlinear  equations  in  one  space  dimension  (to  the  requisite  degree  of  accuracy) 
by  a one-step  noniterative  scheme.  Since  no  iteration  is  required  to  compute  the 
solution  for  a single  time  step,  and  since  only  moderate  effort  is  required  for 
solution  of  the  implicit  difference  equations,  the  method  is  computationally 
efficient;  this  efficiency  is  retained  for  multidimensional  problems  by  using  ADI 
techniques . The  method  is  also  economical  in  terms  of  storage,  in  its  present  form 
requiring  only  two  time-levels  of  storage  for  each  dependent  variable.  Furthermore, 
the  ADI  technique  reduces  multidimensional  problems  to  sequences  of  calculations 
which  are  one-dimensional  in  the  sense  that  easily-solved  narrow  block-banded 
matrices  associated  with  one-dimensional  rows  of  grid  points  are  produced.  Con- 
sequently, only  these  one-dimensional  problems  require  rapid-access  storage  at  any 
given  stage  of  the  solution  procedure,  and  the  remaining  flow  variables  can  be  saved 
on  auxiliary  storage  devices  if  desired. 


4 


J 


R75 -911363-15 


Applicability 

Although  present  attention  is  focused  on  the  compressible  Navier-Stokes 
equations,  the  numerical  method  employed  is  quite  general  and  is  formally  derived 
for  systems  of  governing  equations  which  have  the  following  form: 

dH(4>)dt  = 2H4>)+s(4>)  (i) 

where  0 is  a column  vector  containing  i dependent  variables,  H and  S are  column 
vector  functions  of  0,  and  3 is  a column  vector  whose  elements  are  spatial  differen- 
tial operators  which  may  be  multidimensional.  The  generality  of  Eq.  (l)  allows  the 
method  to  be  developed  concisely  and  permits  various  extensions  and  modifications 
(e.g.,  noncartesian  coordinate  systems,  turbulence  models)  to  be  made  more  or  less 
routinely.  It  should  be  emphasized,  however,  that  the  Jacobian  SH/B0  must  usually 
be  nonsingular  if  the  ADI  techniques  as  applied  to  Eq.  (l)  are  to  be  valid.  A 
necessary  condition  is  that  each  dependent  variable  appear  in  one  or  more  of  the 
governing  equations  as  a time  derivative.  An  exception  would  occur  if  for  instance, 
a variable  having  no  time  derivative  also  appeared  in  only  one  equation,  so  that 
this  equation  could  be  decoupled  from  the  remaining  equations  and  solved  a posteriori 
by  an  alternate  method.  As  a consequence,  the  present  method  is  not  directly  appli- 
cable to  the  incompressible  Navier-Stokes  equations  except  in  one-dimension, 
where  ADI  techniques  are  unnecessary.  For  example,  the  velocity-pressure  form  of 
the  incompressible  equations  has  no  time  derivative  of  pressure,  whereas  the 
vorticity-stream-function  form  has  no  time  derivative  of  stream  function.  For  com- 
puting steady  solutions,  however,  the  addition  of  suitable  "artificial"  time  deri- 
vatives to  the  incompressible  equations,  as  was  done  in  Chorin's  (Ref.  5)  artificial 
compressibility  method,  would  permit  the  application  of  the  present  method.  Alter- 
natively, a low  Mach  number  solution  of  the  compressible  equations  can  be  computed. 
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GOVERNING  EQUATIONS 


The  numerical  method  is  presented  for  flow  in  three  space  dimensions;  two- 
dimensional  problems  can  be  treated  as  a special  case.  For  simplicity,  it  is 
assumed  that  the  fluid  is  a perfect  gas  with  zero  bulk  viscosity  coefficient  and 
constant  molecular  viscosity,  thermal  conductivity,  and  specific  heat.  The 
governing  equations  are  nondimensionalized  by  normalizing  dimensional  variables 
with  the  following  reference  quantities:  distance,  Lr;  velocity,  Ur;  density, 

pr,  temperature;  Tr;  time,  Lr/ur;  enthalpy  Ur2;  and  pressure,  prUr2/g,  where  g is 
the  gravitational  constant.  This  normalization  leads  to  the  following  nondimen- 
sional  parameters:  Mach  number,  M;  Reynolds  number,  Re;  Prandtl  number,  Pr;  and 

specific  heat  ratio,  y.  These  parameters  are  defined  by 


M = ur/C  , Re=/Jr  ur  Lr//JL  , Pr  = cp  /x/k  , X = cp/cv 


(2a-d) 


where  p is  the  molecular  viscosity,  k is  thermal  conductivity,  and  Cp  and  cv  are 
the  specific  heats  at  constant  pressure  and  volume.  The  reference  speed  of  sound, 
c,  is  defined  by  c^  = ygRTr,  where  R is  the  gas  constant. 

With  the  exception  of  the  energy  equation,  the  equations  are  written  in  the 
so-called  conservation  form.  The  foregoing  assumptions  are  convenient  but  not 
essential;  the  treatment  of  alternate  forms  of  the  equations,  arbitrary  equation 
of  state,  and  variable  fluid  properties  is  relatively  straightforward.  With  the 
stated  assumptions,  the  Navier-Stokes  equations  can  be  written  for  Cartesian 
coordinates  (x,y,z)  as  follows:  the  continuity  equation  is 

dp/d t = di-pu)/dx  + d{-p\j)/dy  + di-pw)  /dz  (3a) 

The  momentum  equations  are  represented  by 


d(pTi)/d\  = di-puu)  /<3x  + d(~p  vuO/dy  + di-pwu)  / dz  ~ dp/d'x  + F (3h) 


The  energy  equation  is 
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In  Eqs . (3).,  u,  v,  and  w are  the  x,  y and  z components  of  the  velocity  vector,  U; 
p is  density;  T is  temperature;  p is  pressure;  t is  time,  and  v is  the  gradient 
operator.  The  symbols  u,  x denote  u,  x;  v,  y;  w,  z,  respectively,  for  the  x,  y, 
and  z momentum  equations. 


The  force  due  to  viscous  stress,  F,  is  given  by 


F = 


R e 


[a2u/ax2  + a2u  /ay2  +a2u/az2  + -y-  a(v-u) /ax 


(4) 


The  dissipation  function,  ?,  is  given  by 


0 = 2 


au 

dx 


av  \ . / aw  \ ] . / au  av  \ 
-TTl  +[-d7)  J ('ay  + dx  ) 


+ +(4sl+  Js.\* _ 2 n . 


az 


ay 


az 


ax  / 


v-u 


)' 


(5) 


The  pressure  can  be  eliminated  as  a dependent  variable  by  means  of  the  equation 
of  state  for  a perfect  gas, 


p=pT/YMZ  (6) 

The  continuity,  momentum  and  energy  equations  thus  constitute  a system  of  five 
equations  in  the  dependent  variables  p,  u,  v,  w,  and  T.  The  definition  of  total 
enthalpy  E is 

E=T/(r-l)M2  +q2/2  (7) 

o p 2 p 

where  q11-  = u + v + w . In  numerous  problems  of  interest,  it  can  be  assumed  that 
the  total  enthalpy  is  a constant  EQ  provided  there  is  no  heat  addition.  This 
assumption  is  reasonable  for  inviscid  flow  regions  with  or  without  shocks  and  for 
boundary  layers  on  adiabatic  walls  provided  the  Prandtl  number  is  unity.  In  this 
circumstance,  Eqs.  (6-7)  can  be  combined  to  produce  an  adiabatic  equation  of  state, 

p =/°(E  0 -q2/2)  (x-l)  /X  (8) 

If  pressure  is  eliminated  in  the  momentum  equations  by  means  of  (8),  then  solution 
of  the  energy  equation  (3c)  is  unnecessary,  and  a significant  reduction  in  compu- 
tational effort  is  effected.  The  temperature  field  is  then  determined  a posteriori 
from  Eq.  (7).  This  simplification,  although  convenient  and  available  as  an  option, 
was  not  used  for  any  of  the  calculations  presented  here  for  laminar  flow. 
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NUMERICAL  METHOD 


Previous  Work 

Although  several  methods  based  on  implicit  schemes  have  been  developed  for 
incompressible  flows  (e.g.,  Pearson,  Ref.  6;  Chorin,  Ref.  1) , most  previous  methods 
for  the  compressible  Navier-Stokes  equations  have  employed  explicit  schemes. 
Nevertheless,  a semi-implicit  method  has  been  developed  by  Harlow  and  Amsden  (Ref. 

8)  for  use  over  the  entire  spectrum  of  Mach  numbers  from  incompressible  to  hyper- 
sonic. However,  the  Harlow-Amsden  method  treats  viscous  terms  explicitly  and, 
unlike  alternating-direction  methods,  requires  the  solution  of  multidimensional 
implicit  difference  equations,  which  tends  to  be  time  consuming.  In  an  independent 
investigation,  Baum  and  Ndefo  (Ref.  ^)  developed  a two-dimensional  implicit  method 
which  is  patterned  after  the  original  Peaceman-Rachford  (Ref.  10)  ADI  technique. 
Perhaps  the  most  significant  difference  between  the  Baum-Ndefo  and  present 
methods  is  that  the  Baum-Ndefo  method  employs  iterative  techniques,  solving  non- 
linear difference  equations  as  a sequence  of  linear  equation,  whereas  in  the  pre- 
sent method,  the  difference  equations  are  linearized  about  the  solution  at  the  pre- 
vious time  step  and  solved  without  iteration.  In  principle,  the  solution  of  non- 
linear difference  equations  is  of  course  attractive,  as  this  removes  any  limitations 
which  might  arise  from  the  linearization  process.  It  is,  however,  a time-consuming 
process  to  solve  the  nonlinear  difference  equations,  and  as  far  as  temporal  accuracy 
is  concerned,  the  additional  computational  effort  required  by  the  solution  of 
nonlinear  difference  equations  might  be  as  well  spent  by  reducing  the  time  step 
and  proceeding  with  a satisfactory  linearization  (McDonald  & Briley,  Ref.  11). 

On  the  other  hand,  if  a steady  solution  is  the  only  objective,  then  temporal 
accuracy  is  of  little  concern,  and  a stable  method  requiring  minimal  computational 
effort  per  time  step  (the  present  objective)  is  attractive.  The  topic  of  nonlinear 
truncation  errors  is  discussed  further  by  McDonald  and  Briley  (Ref.  11). 

The  present  numerical  method  was  developed  for  the  Navier-Stokes  equations 
by  Briley  and  McDonald  (Ref.  12)  and  for  steady  supersonic  flows  by  McDonald  and 
Briley  (Ref.  11).  Here  the  method  is  formalized  for  mixed  parabolic-hyperbolic 
systems  having  the  form  of  Eq.  (l),  and  is  applied  to  the  Navier-Stokes  equations. 

Recently,  Beam  and  Warming  (Ref.  13)  have  specialized  this  same  linearization 
procedure  to  equations  in  conservation-law  form,  with  emphasis  on  the  inviscid 
Euler  equations.  Beam  and  Warming  also  employed  compact  spatial  differencing 
techniques  analogous  to  those  discussed  by  Mitchell  (Ref.  l4,  p.  51)  but  in  more 
general  circumstances.  Although  highly  effective  as  applied  by  Beam  and  Warming  to 
conservation  laws,  the  compact  spatial  differencing  appears  to  lose  its  advantage 
in  computational  efficiency  over  more  conventional  differencing  techniques,  when- 
ever both  first  and  second-order  derivatives  are  present.  This  point  will  be 
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discussed  subsequently.  Beam  and  Warming  also  suggested  valuable  techniques  for 
use  in  computing  flows  containing  shocks. 


Linearization  Technique 


Background 


A number  of  techniques  have  been  used  for  implicit  solution  of  the  following 
first-order  nonlinear  scalar  equation  in  one  dependent  variable  0(x,t): 

d<f>/d\  = dG(<£)/<3x  (9) 

Special  cases  of  Eq.  (9)  include  the  conservation  form  if  F(0)  = 1,  and  quasilinear 
form  if  G(0)  = 0.  Previous  implicit  methods  for  Eq.  (9)  which  employ  nonlinear 
difference  equations  and  also  methods  based  on  two-step  predictor-corrector  schemes 
are  discussed  by  Ames  (Ref.  15,  p.  82)  and  von  Rosenberg  (Ref.  l6,  p.  56).  One 
such  method  is  to  difference  nonlinear  nonlinear  terms  directly  at  the  implicit 
time  level  to  obtain  nonlinear  implicit  difference  equations;  these  are  then  solved 
iteratively  by  a procedure  such  as  Newton's  method.  Although  otherwise  attractive, 
there  may  be  difficulty  with  convergence  in  the  iterative  solution  of  the  nonlinear 
difference  equations,  and  some  efficiency  is  sacrificed  by  the  need  for  iteration. 

An  implicit  predictor-corrector  technique  has  been  devised  by  Douglas  and  Jones 
(Ref.  17)  which  is  applicable  to  the  quasilinear  case  (G  = 0)  of  Eq.  (9).  The 
first  step  of  their  procedure  is  to  linearize  the  equation  by  evaluating  the  non- 
linear coefficient  as  F(0n)  and  to  predict  values  of  0n+' 2 using  either  the  backward 
difference  or  the  Crank-Nijolson  scheme.  Values  for  0n+1  are  then  computed  in  a 
similar  manner  using  F(0n+2')  and  the  Crahk-Nicolson  scheme.  Gourlay  and  Morris 
(Ref.  18)  have  also  proposed  implicit  predictor- corrector  techniques  which  can  be 
applied  to  Eq.  (9).  In  the  conservative  case  (F=l),  their  technique  is  to  define 
G(0)  by  the  relation  G(0)  = 0G(0)  when  such  a definition  exists,  and  to  evaluate 
G(0n+1)  -using  values  for  0n+I  computed  by  an  explicit  predictor  scheme.  With  G 
thereby  known  at  the  implicit  time  level,  the  equation  can  be  treated  as  linear, 
and  corrected  values  of  0n+-'-  are  computed  by  the  Crank-Nicolson  scheme. 

A technique  is  described  here  for  deriving  linear  implicit  difference 
approximations  for  nonlinear  differential  equations.  The  technique  is  based  on  an 
expansion  of  nonlinear  implicit  terms  about  the  solution  at  the  known  time  level, 
tn,  and  leads  to  a one-step,  two-level  scheme  which,  being  linear  in  unknown 
(implicit)  quantities,  can  be  solved  efficiently  without  iteration.  This  idea  was 
applied  by  Richtmyer  and  Morton  (Ref.  19,  p.  203)  to  a scalar  nonlinear  diffusion 
equation.  Here,  the  technique  is  developed  for  problems  governed  by  & nonlinear 
equations  in  & dependent  variables  which  are  functions  of  time  and  space  coordinates. 
Attention  is  restricted  to  nonlinear  systems  having  the  form  of  Eq.  (l). 
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The  solution  domain  is  discretized  by  grid  points  having  equal^spacings , 

Ax,  Ay,  and  Az,  in  the  x,  y,  and  z directions,  respectively,  and  an  arbitrary  time 
step,  At.  Provisions  for  nonuniform  grid  spacing  will  be  introduced  subsequently. 
The  subscripts  i,  j,  k and  superscript  n are  grid  point  indices  associated  with  x, 
y,  z,  and  t,  respectively,  and  thus  ^ denotes  } (xq,yj  . It  is  assumed 

that  the  solution  is  known  at  the  n level,  tn,  and  is  desired  at  the  (n+l)  level, 
tn+l.  At  the  risk  of  an  occasional  ambiguity,  one  or  more  of  the  subscripts  is  fre- 
quently omitted,  so  that  0n  is  equivalent  to  0?  . , . 

Linearized,  Difference  Scheme 

The  linearized  difference  approximation  is  derived  from  the  following  implicit 
time-difference  replacement  of  Eq.  (l): 


(Hn  + '-Hn)/At=/3 


'2>  (<£n+l)  + Sn  + l]+(i-/3) 


, n , n 
i<p  ) S 


(10) 


where,  for  example,  Hn+‘*~  = H(0n+'*').  The  form  of  D and  the  spatial  differencing 
are  as  yet  unspecified.  A parameter  p(0  < 3 < l)  has  been  introduced  so  as  to 
permit  a variable  centering  of  the  scheme  in  time.  Equation  (10)  produces  a backward 
difference  formulation  for  p = 1 and  a Crank-Nicolson  formulation  for  p - l/2. 
Unconditional  stability  is  anticipated  for  p > l/2. 

The  linearization  is  performed  by  a two-step  process  of  expansion  about  the 
known  time  level  tn  and  subsequent  approximation  of  the  quantity  (S0/dt)nAt,  which 
arises  from  chain  rule  differentiation,  by  (0n+^--0n).  The  result  is 

Hn+I  = Hn  +(dH/d<£)n  {4>n  + '-cf>n)  + 0(At)2 

Sn+l=Sn  + Us/#)n  (<£n  + l -<£  n)  +0  ( At)2 

■2>(^>n+l)  = 2)(<f>°)+{d  ^/<3<^)(<^n+,-</>n)  + 0(At)2 

The  matrices  3h/50  and  dS / 90  are  standard  Jacobians  whose  elements  are  defined, 
for  example,  by  (gH/50)  = 3Hq/S0r.  The  operator  elements  of  the  matrix  S2>/ 30 

are  similarly  ordered,  i.e.,  (32>/S0)qr  - d-^q/^r’  however,  the  intended  meaning 
of  the  operator  elements  requires  some  clarification.  For  the  qth  row,  the  opera- 
tion (d-2q/c$)n  (0n+1  - 0n)  is  understood  to  mean  that  { s/ dt  [0(x,y,z,t)]  }n  At 
is  computed  and  that  all  occurrences  of  (30r/3t)n  arising  from  chain  rule  differen- 
tiation are  replaced  by  (0^+^  - 0r)/At. 

After  linearization  as  in  Eqs.  (ll),  Eq.  (10)  becomes  the  following  linear 
implicit  time-differenced  scheme: 

(dHn/a<£)(<£n  + l-<£n) /At  =2)(4>n)  +sn  + /3  {32)  /d<f>  + dsn /d4>){<f>n+'-4>n)  (12) 


(11a) 

(lib) 

(11c) 
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Although  Hn+^~  is  linearized  to  second  order  in  Eq.  (lla),  the  division  Toy  At  in 
Eq.  (10)  introduces  an  error  term  of  order  At.  A technique  for  maintaining  formal 
second-order  accuracy  in  the  presence  of  nonlinear  time  derivatives  is  discussed 
by  McDonald  and  Briley  (Ref.  11),  however  a three-level  scheme  results.  Second- 
order  temporal  accuracy  can  also  be  obtained  (for  3 = l/2)  by  a change  in  dependent 
variable  to  0 = H (0),  provided  this  is  convenient,  since  the  nonlinear  time  deriva- 
tive is  then  eliminated.  The  temporal  accuracy  is  independent  of  the  spatial  accuracy. 


On  examination,  it  can  be  seen  that  Eq.  (12)  is  linear  in  the  quantity 
(0n+l  _ 0n)  an^  that  all  other  quantities  are  either  known  or  evaluated  at  the  n 
level.  Computationally,  it  is  convenient  to  solve  Eq.  (12)  for  (0n+l  - 0n)  rather 
than  0n+^.  This  both  simplifies  Eq.  (12)  and  reduces  roundoff  errors,  since  it  is 
presumably  better  to  compute  a small  O(At)  change  in  an  0(l)  quantity  than  the 
quantity  itself.  To  simplify  the  notation,  a new  dependent  variable  i|f  defined  by 


'P  = 4>-4>n  (13) 

is  introduced,  and  thus  = 0n+1  - 0n,  and  \|rn  = 0.  It  is  also  convenient  to 

rewrite  Eq.  (12)  in  the  following  simplified  form: 


(A+ At^)^n  + ' 


At 


JM4>n)+s" 


(l4a) 


where  the  following  symbols  have  been  introduced  to  simplify  the  notation: 

As  dHn/di> -/3At(dSn/'<3<£)  < 

J=-P(d  %/d<fi) 

It  is  noted  that  i(i|/)  is  a linear  transformation  and  thus  f(0)  = 0.  Furthermore, 
if  -2>  (0)  is  linear,  then  f(i|;)  = -(32>(i|r). 


(14b) 

(14c) 


Spatial  differencing  of  Eq.  (l4a)  is  accomplished  simply  by  replacing  derivative 
operators  such  as  d/dx,  S2/ by  corresponding  finite  difference  operators,  Dx,  Dx. 
Henceforth,  it  is  assumed  that  2>  and  £ have  been  discretized  in  this  manner,  unless 
otherwise  noted. 


In  practice,  difference  schemes  generated  by  the  foregoing  procedures  are  often 
quite  simple.  For  example,  the  continuity  equation  (3a)  becomes 


(pn  + '-pn)  = -P  At 


I D~f/+lunV“n  + ,+(^)/unl 

(u,x)  X L P J 


, . , n n+l  ~n+l 
which  couples  p and  u 


(15) 
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Discussion 


Before  proceeding,  some  general  observations  seem  appropriate.  The  foregoing 
linearization  technique  assumes  only  Taylor  expandability,  an  assumption  already 
implicit  in  the  use  of  a finite  difference  method.  The  governing  equations  and 
boundary  conditions  are  addressed  directly  as  a system  of  coupled  nonlinear  equations 
which  collectively  determine  the  solution.  The  approach  thus  seems  more  natural 
than  that  of  making  ad  hoc  linearization  and  decoupling  approximations,  as  is  often 
done  in  applying  implicit  schemes  to  coupled  and/or  nonlinear  partial  differential 
equations.  With  the  present  approach,  it  is  not  necessary  to  associate  each 
governing  equation  and  boundary  condition  with  a particular  dependent  variable 
(e.g.,  assumed  u is  governed  by  the  x momentum  equation,  p by  continuity,  etc.) 
and  then  to  identify  various  "nonlinear  coefficients"  and  "coupling  terms"  which 
must  then  be  treated  by  lagging,  predictor- corrector  techniques,  or  iteration.  The 
Taylor  expansion  procedure  is  analogous  to  that  used  in  the  generalized  Newton- 
Raphson  or  quasilinearization  methods  for  iterative  solution  of  nonlinear  systems 
by  expansion  about  a known  current  guess  at  the  solution  (e.g.,  Bellman  & Kalaba, 
Ref.  20).  However,  the  concept  of  expanding  about  the  previous  time  level  has 
apparently  not  been  employed  to  produce  a noniterative  implicit  time-dependent 
scheme  for  coupled  equations,  wherein  nonlinear  terms  are  approximated  to  a level 
of  accuracy  commensurate  with  that  of  the  time  differencing.  The  linearization 
technique  also  permits  the  implicit  treatment  of  coupled  nonlinear  boundary  condi- 
tions, such  as  stagnation  pressure  and  enthalpy  at  subsonic  inlet  boundaries,  and 
in  practice,  this  latter  feature  was  found  to  be  crucial  to  the  stability  of  the 
overall  method. 


Application  of  Alternating-Direction  Techniques 

Solution  of  Eq.  (l4a)  is  accomplished  by  application  of  an  alternating- 
direction  implicit  (ADl)  technique  for  parabolic-hyperbolic  equations.  The  original 
ADI  method  was  introduced  by  Peaceman  and  Rachford  (Ref.  10)  and  Douglas  (Ref.  21); 
however,  the  alternating-direction  concept  has  since  been  expanded  and  generalized. 

A discussion  of  various  alternating-direction  techniques  is  given  by  Mitchell  (Ref. 
l4)  and  Yanenko  (Ref.  22). 

The  present  technique  is  simply  an  application  of  the  very  general  procedure 
developed  by  Douglas  and  Gunn  (Ref.  4)  for  generating  ADI  schemes  as  perturbations 
of  fundamental  implicit  difference  schemes  such  as  the  backward-difference  or 
Crahk-Nicolson  schemes. 

For  the  present,  it  will  be  assumed  that  2>(0)  contains  derivatives  of  first 
and  second  order  with  respect  to  x,  y,  and  z,  but  no  mixed  derivatives.  In  this 
case,  2>  can  be  split  into  three  operators,  Y>x,  2>y,  2>z  associated  with  the  x,  y,  and 
z coordinates  and  each  having  the  functional  formY>x  - *$?(0,  b/ Sx,  B2/&>?)  for  a 
typical  coordinate  x.  Equation  (l4a)  then  becomes 
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_A+At(/X+/  y+/Z) 


'/'n  + l = At[(J2> 


x+3>y+^z)<£n+Sn] 


(16) 


Recalling  that  _/(\jtn)  = 0,  the  Douglas -Gunn  representation  of  Eq.  (l6)  can  be 
written  as  the  following  three-step  solution  procedure: 


(A+At/x)^:,t  = At[(2)x+^y  + ^z)^n+Sn] 

(17a) 

(A  + At  JyW**  = At* 

(17b) 

(A  + At  Jz)fn  + i 

(17c) 

where  i|r*  and  ij/**  are  intermediate  solutions.  It  will  be  shown  subsequently  that 
each  of  Eqs.  (17)  can  be  written  in  narrow  block-banded  matrix  form  and  solved 
by  efficient  block- elimination  methods.  If  and  \|f**  are  eliminated,  Eqs.  (17) 
become 


(A  + At  JK) 


(A  + At  if  y ) A" 


(A+At  I z)^ 


n + 1 


= At 


^y  + 


3>z> 


4>"+sn] 


(18) 


If  the  multiplication  on  the  left-hand  side  of  Eq.  (18)  is  performed,  it  becomes 
apparent  that  Eq.  (l8)  approximates  Eq.  (l6)  to  order  (At)2.  Although  the  stability 
of  Eqs.  (17)  has  not  been  established  in  circumstances  sufficiently  general  to 
encompass  the  Navier-Stokes  equations,  it  is  often  suggested  (e.g.,  Richtmyer  & 
Morton,  Ref.  19,  p.  215)  that  the  scheme  is  stable  and  accurate  under  conditions 
more  general  than  those  for  which  rigorous  proofs  are  available.  This  latter  notion 
is  adopted  here  as  a working  hypothesis  supported  by  favorable  results  obtained 
in  actual  computations . 

A major  attraction  of  the  Douglas -Gunn  scheme  is  that  the  intermediate 
solutions  i|f*  and  ijr  are  consistent  approximations  to  i|rn+l.  Furthermore,  for  steady 
solutions,  i|rn  = ijr*  = i|/**  = \|/n+2-  independent  of  At.  Thus,  physical  boundary  condi- 
tions for  ij/n+1  can  be  used  in  the  intermediate  steps  without  a serious  loss  in 
accuracy  and  with  no  loss  for  steady  solutions.  In  this  respect,  the  Douglas -Gunn 
scheme  appears  to  have  an  advantage  over  locally  one- dimensional  (LOD)  or  "splitting" 
schemes,  and  other  schemes  whose  intermediate  steps  do  not  satisfy  the  consistency 
condition.  The  lack  of  consistency  in  the  intermediate  steps  complicates  the  treat- 
ment of  boundary  conditions  and,  according  to  Yanenko  (Ref.  22,  p.  33) 5 does  not 
permit  the  use  of  asymptotically  large  time  steps.  It  is  not  clear  that  this  advan- 
tage of  the  Douglas-Gunn  scheme  would  always  outweigh  other  benefits  which  might 
be  derived  from  an  alternative  scheme.  However,  since  the  ADI  scheme  can  be  viewed 
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as  an  approximate  technique  for  solving  the  fundamental  difference  scheme,  Eq. 
(l4a),  alternate  ADI  schemes  can  readily  be  used  within  the  present  formulation. 

It  is  worth  noting  that  the  operator  2>  can  be  split  into  any  number  of 
components  which  need  not  be  associated  with  a particular  coordinate  direction. 

As  pointed  out  by  Douglas  and  Gunn  (Ref.  4),  the  criterion  for  identifying  sub- 
operators is  that  the  associated  matrices  be  "easily  solved"  (i.e.  narrow-banded). 
Thus,  mixed  derivatives  can  be  treated  implicitly  within  the  ADI  framework, 
although  this  would  increase  the  number  of  intermediate  steps  and  thereby  compli- 
cate the  solution  procedure.  Finally,  only  minor  changes  are  introduced  if,  in 
the  foregoing  development  of  the  numerical  method,  H,  2>,  and  S are  functions  of 
the  spatial  coordinates  and  time,  as  well  as  0. 
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Solution  of  the  Implicit  Difference  Equations 


Second-Order  Spatial  Differences 


Since  each  of  Eqs.  (17)  is  implicit  in  only  one  coordinate  direction,  the 
solution  procedure  can  be  discussed  with  reference  to  a one -dimensional  problem. 

For  simplicity,  it  is  sufficient  to  consider  Eq.  (17a)  with  2>  ,2>z  = 0.  For  the 
moment,  attention  is  focused  on  the  following  three-point  difference  formulas: 

=[aA_  + (l-  a)  A+]<£/Ax  = (<3<£/d'5?)j  + 0 [a3(  2+(a -1/2)  Ax]  (19a) 

d/<£MA+A_)<MA'x)2  = (a2^/<Jx2)i+0(A'5if2)  (19b) 


for  a typical  coordinate  y.  Here,  A_=  fa  - fa_i>  A+  = 0i+1  - <t>±>  and  a parameter 
ot  has  been  introduced  (0  <;  a £ l)  so  as  to  permit  continuous  variation  from  back- 
ward to  forward  differences.  The  standard  central  difference  formula  is  recovered 
for  a = \ and  was  used  for  all  solutions  reported  here. 

As  an  example,  suppose  that  the  qth  component  of  3>x  has  the  form 

3 A ^ 

"ax"  Glq^)  + ^2q  1 GZq($>)  (20) 


where  F and  G are  column  vector  functions  having  the  same  but  an  arbitrary  number 
of  components;  F^  denotes  the  transpose  of  F.  The  form  of  Eq.  (20)  permits  govern- 
ing equations  having  any  number  of  first  and  second  derivative  terms.  Then, 


(a3xq/w<*"+'-.»n)sFT_l 


dG 


iq 


(qbn+l-<*>n)  + 


ilk 

d<f> 


<+n+l-*n) 


dG 


iq 


<3  x 


+ F, 


d2  dG2q 

2q  dx2  ~d^~ 


dF, 


(*n+l-*n)^ 


(<pn+l-<pn) 


2 

dG 


(21) 


dx 


2q 

~T~ 


It  is  now  possible  to  describe  the  solution  procedure  for  Eq.  (l7a)  for  the 
one -dimensional  case  with  2>v  given  by  Eq.  (20)  and  difference  formulas  given  by 
Eq,  (19).  Because  of  the  spatial  difference  operators,  Dx  and  Dx  , Eq.  (17a) 
contains  \|/t  tyt,  and  consequently,  the  system  of  linear  equations  generated 

by  writing  Eq.  (17a)  at  successive  grid  points  x^  can  be  written  in  block- tridiagonal 
form  (simple  tridiagonal  for  scalar  equations,  i = l).  The  block-tridiagonal 
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matrix  structure  emerges  from  rewriting  Eq.  (17a)  as 


an  i/t*  + +c?f* 

i i - 1 ti  *i+i 


(22) 


where  a,  b,  c are  square  matrices  and  d is  a column  vector,  each  containing  only 
n-level  quantities.  The  qth  component  of  d and  qth  row-components  of  a,  b,  c are 
given  by 


(<>!Y  [WF|Tq/a+)"b|q).n_  + (F^)"(2GIq/2(#,)7|  ]a/Ax 

-[(dFT/a+lf  (o2,)"_,  +(FaTq)^(aGs,q/d<#>)"_1]  /(A*)’ 


(23a) 


( bn j )q  = (dHq/d<t>)-  //3At  - (asq/a<£)" 

+ KV*K<G|q)"  +(F,F)"  <dG,q /a<#*)"]  (l-2a)/Ax 
+a[(dFjq/a*)J(G2q),  +(F2Tq)"(aG2q/a*)"]/(Ax)2 


(cpq  = [(aF,Tq/a^)|(G,q)"+|  +(F,q)|  (aGlq/a<#>)"+l]  (a-o/Ax  (230) 

- [(ap|q  /a*)1)  (s2q), +l  +(F2Tql"  (aG2q/a</>)"+l]  /( Ax)2 

r , (23d) 

(d")a=  [,S(/l)  + sn]//3 


When  applied  at  successive  grid  points,  Eq.  (22)  generates  a block-tridiagonal  system 
of  equations  for  i|t*  which,  after  appropriate  treatment  of  boundary  conditions,  can 
be  solved  efficiently  using  standard  block-elimination  methods  as  discussed  by 
Isaacson  and  Keller  (Ref.  23,  p.  58).  The  solution  procedure  for  Eqs.  (I7b&c)  is 
analogous  to  that  just  described  for  Eq.  (17a).  It  is  worth  noting  that  the  spatial 
difference  parameter  <y  canbe  varied  with  i or  even  term  by  term.  For  example,  an 
"upwind  difference"  formula  can  be  obtained  if  a is  chosen  as  1 or  -1  depending  on 
the  sign  of  the  elements  of  Fj_. 
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Fourth-Order  Spatial  Differences 


Fourth-order  spatial  accuracy  can  be  obtained  by  using  the  following  standard 
five-point  difference  formulas  in  place  of  Eqs.  (19) : 

D^4>  = IhAir  ~6~A+  A_KA++A_)  <t>  =(d^/ax)i+o(Ax4)  (24a) 

d£</>  = — {—  (i--rk-  a ,a_)(a,  a )<t>  =(a2^>/ax2)i+o(Ax4)  (24b) 

(Ax) 

In  this  instance,  the  block-tridiagonal  structure  of  Eq.  (22)  expands  to  block  quin- 
diagonal.  Block-quindiagonal  systems  are  easily  solved  using  banded  Gaussian  block- 
elimination.  The  question  arises,  however,  whether  the  increased  accuracy  thus 
obtained  is  worth  the  additional  computation  involved  in  solving  block-quindiagonal 
systems.  This  topic  will  be  discussed  subsequently. 

Recently,  there  has  been  revived  interest  (Orszag  & Israeli,  Ref.  l)  in  the 
following  fourth-order  compact  difference  formulas  involving  only  three  grid  points : 


= 


iw?  irfH-r* = <«*/«>,  +o(a'x4) 

2Ax  (|  + -£-A  A_)  1 


D*<fi  = -L-Z  {A+A--\ £ ={dz4>/frjfr  + 0(Ax4) 

(Axf  (l  + ^-A+A_)  1 


(25a) 

(25b) 


To  implement  Eqs.  (25)  within  the  ADI  framework,  the  difference  equations  are  multi- 
plied by  the  denominator  of  Eq.  (25a)  or  (25b)  at  the  appropriate  step  in  the  solu- 
tion procedure  (cf.  Mitchell,  Ref.  l4;  Beam  & Warming,  Ref.  13). 

The  attraction  of  Eqs.  (25)  over  (24)  is  that  Eqs.  (25)  lead  to  block-tridiagonal 
rather  than  block-quindiagonal  equations.  One  difficulty  which  arises  with  Eqs.  (25), 
however,  is  that  the  denominators  are  different  and  must  be  removed  separately  when 
both  first  and  second  order  derivatives  are  present,  as  in  the  Navier-Stokes 
equations.  Consequently,  first  and  second  derivative  terms  require  a separate  step 
in  the  ADI  procedure,  and  thus  two  block-tridiagonal  systems  must  be  solved  for 
Eqs.  (25),  as  opposed  to  one  block-quindiagonal  system  for  Eqs.  (24),  in  each 
coordinate  direction  that  both  first  and  second  derivatives  appear.  In  terms  of 
computational  effort,  these  requirements  are  approximately  equal,  as  will  be  shown 
in  the  following  section,  and  from  this  standpoint,  there  is  little  to  choose  between 
Eqs.  (24)  & (25).  Both  Eqs.  (24)  and  (25)  suffer  from  complications  in  applying 
boundary  conditions,  as  is  usually  the  case  when  a difference  formula  is  of  higher 
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order  than  the  derivative  approximated.  One  disadvantage  of  Eqs . (25)  is  that  the 
evaluation  of  explicit  derivatives  requires  the  solution  of  a simple  tridiagonal 
system  for  each  dependent  variable  present.  Another  limitation  of  Eqs.  (25)  arises 
when  terms  containing  derivatives  with  nonconstant  coefficients  are  present,  for 
example,  variable  viscosity  terms.  Each  nonconstant  coefficient  can  be  removed  by 
division  prior  to  clearing  the  denominator  of  Eqs.  (25),  but  as  a consequence,  each 
term  of  this  type  in  the  same  governing  equation  also  requires  a separate  step  in 
the  ADI  procedure,  and  the  solution  of  a block-tridiagonal  system.  The  authors  are 
thus  inclined  to  favor  the  five-point  formulas  Eqs.  (24),  since  their  use  and  effi- 
ciency is  less  dependent  on  the  form  of  the  governing  equations. 


Computing  Requirements 

Various  block-elimination  algorithms  can  be  devised  for  solution  of  equations 
with  block-banded  matrix  structures  (cf.,  Isaacson  & Keller,  Ref.  23).  Such  algo- 
rithms can  be  derived  using  variants  of  Gaussian  elimination  for  a banded  matrix, 
but  with  the  square  submatrix  elements  of  the  banded  matrix  processed  using  matrix 
algebra.  Thus,  operations  involving  matrix  subelements  are  not  assumed  to  commute, 
and  division  by  a matrix  subelement  is  accomplished  by  computing  the  inverse  and 
multiplying.  Following  this  procedure,  the  authors  have  developed  algorithms  for 
both  block-tridiagonal  and  block-quindiagonal  systems  arising  from  the  respective 
second  and  fourth-order  difference  formulas,  Eqs.  (19)  & (24).  Each  algorithm 
requires  only  one  inverse  per  grid  point.  A standard  operation  count  (scalar 
multiplications  and  divisions)  has  been  performed  for  systems  with  L x L block 
elements  and  N diagonal  block  elements,  i.e.,  L coupled  equations  along  N grid 
points.  The  block- tridiagonal  scheme  requires  (3N-2)  (L^  + L2)  operations,  the  same 
as  the  matrix  factorization  scheme  of  Isaacson  and  KeILer(Ref.  23);  the  block-quin- 
diagonal scheme  requires  (7N-10)  l3  +(5N-6)l£  operations,  which  is  only  slightly 
more  than  twice  the  block-tridiagonal  number. 

Orszag  and  Israeli  (Ref.  1,  p.  284)  have  estimated  that  fourth-order  schemes 
achieve  results  in  the  5 percent  accuracy  range  with  approximately  half  the  number 
of  grid  points  in  each  direction,  as  compared  with  second-order  schemes.  Thus  based 
on  the  foregoing  estimate  and  operation  counts,. it  is  concluded  that  for  one-dimensional 
problems,  use  of  the  fourth-order  scheme  is  roughly  equivalent  in  terms  of  accuracy 
and  computational  cost  to  using  the  second-order  scheme  with  twice  as  many  grid  points. 
In  two  and  three  dimensions ,. however , the  second-order  scheme  requires  four  and  eight 
times  as  many  grid  points,  respectively,  to  obtain  accuracy  comparable  to  fourth- 
order  schemes,  and  the  fourth -order  scheme  is  well  worth  the  factor  of  two  in  com- 
putational effort  per  grid  point. 

Assuming  there  are  E grid  points  in  each  coordinate  direction,  the  total  num- 
ber of  operations  for  a single  time  step  is  obtained  from  the  operation  count  for 
solution  of  one  block-banded  system  by  multiplying  by  2N  and  3N2  for  two  and  three 
dimensions,  respectively. 
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For  the  particular  case  of  the  Navier-Stokes  equations  (3)  with  p eliminated 
using  Eq.  (6),  it  is  possible  to  reduce  the  computational  effort  substantially  by 
taking  advantage  of  the  special  nature  of  the  coupling  during  each  ADI  sweep.  In 
this  case,  it  is  only  necessary  to  solve  one  block-banded  system  with  L = 3,  as 
well  as  two  simple  banded  systems  (L  = l).  This  can  be  seen  by  careful  examination 
of  Eqs . (3).  During  the  first  step  of  the  ADI  procedure,  only  derivatives  with 
respect  to  x and  t from  Eqs.  (3)  appear  in  the  implicit  difference  equations.  In 
the  continuity,  x momentum,  and  energy  equations,  these  implicitly  treated  terms 
contain  p,  u,  and  T,  but  not  v and  w;  therefore,  the  difference  equations  from  these 
three  equations  can  be  solved  for  p*,  u*,  and  T*  during  the  first  ADI  step  by  sol- 
ving a block -banded  system  with  L = 3.  Having  obtained  values  for  p*,  u*,  and  T* 
in  this  manner,  the  difference  equations  for  the  y and  z momentum  equations  can  then 
be  solved  independently  for  v*  and  w*;  since  the  y and  z momentum  equations  are 
uncoupled  with  respect  to  v*  and  w*  during  this  ADI  step,  the  latter  computation 
only  requires  the  solution  of  two  simple  banded  systems  (L  = l).  A similar  situation 
exists  for  the  remaining  two  steps  of  the  ADI  procedure,  except  that  during  the 
second  step  (y  derivatives  treated  implicitly),  the  difference  equations  from  the 
continuity,  y momentum,  and  energy  equations  are  solved  as  coupled  equations  for 
p**,  v** , and  T**,  and  during  the  third  step  (z  derivatives  treated  implicitly),  the 
difference  equations  from  the  continuity,  z momentum,  and  energy  equations  are  solved 
as  coupled  equations  for  pn+^,  w114^,  and  Tn+-'-. 

For  tridiagonal  systems,  the  operation  count  is  reduced  from  order  450  N for 
one  L = 5 system  to  order  118  N for  one  L = 3 and  two  L = 1 systems.  For  quindiag- 
onal  systems,  these  estimates  are  1000  H and  258  N,  respectively.  Consequently,  the 
arrangement  leading  to  three  coupled  and  two  uncoupled  equations  is  quite  worthwhile. 
For  comparison,  it  is  noted  that  in  the  case  of  the  Navier-Stokes  equations  (3), 
merely  evaluating  the  right-hand  side  of  Eq.  (17a),  which  would  be  a minimum  require- 
ment for  a one-step  explicit  scheme,  requires  302  N operations  for  a 3-point  differ- 
ence formula  and  488  N operations  for  a 5-P°int  formula. 

In  view  of  the  many  factors  involved,  it  is  difficult  to  evaluate  precisely 
or  with  any  generality  the  overall  computational  efficiency  of  the  present  method 
relative  to  various  other  methods.  However,  the  foregoing  operational  counts  show 
that  the  effort  expended  to  solve  the  implicit  difference  equations  by  block-elimi- 
nation is  not  excessive  compared  with  that  necessary  simply  to  evaluate  the  differ- 
enced Navier-Stokes  equations,  let  alone  the  various  other  bookkeeping  tasks  present 
in  most  large-scale  computer  programs  for  fluid  dynamics  problems.  In  the  solutions 
presented  here,  the  solution  of  the  tridiagonal  and  block-tridiagonal  systems  using 
double  precision  arithmetic  required  only  about  one  third  to  one  half  of  the  total 
compirter  time  per  time  step. 
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APPLICATION  TO  FLOW  IN  A STRAIGHT  DUCT 


Problem  Formulation  and  Boundary  Conditions 


To  explore  the  stability  properties  and  general  capabilities  of  the  present 
method,  some  test  calculations  were  made  for  three-dimensional  laminar  subsonic 
flow  in  the  entrance  region  of  a straight  duct  with  rectangular  cross  section. 

The  straight  duct  geometry  was  chosen  for  its  simplicity,  and  consideration  was 
limited  to  subsonic  flow  to  avoid  the  possible  occurrence  of  shocks.  A consider- 
ation of  shocks  and  other  complicating  factors  which  would  have  hindered  the 
orderly  development  of  the  method  is  reserved  for  future  studies.  Application 
of  the  present  method  is  straightforward  once  H,  S,  and  2>  are  identified,  and 
for  the  present  calculations,  these  quantities  are  given  in  the  Appendix.  In 
applying  the  numerical  method,  the  dissipation  term,  $,  defined  by  Eq.  (5),  and 
the  viscous  terms  in  Eq.  (4)  containing  v-U  were  treated  explicitly  by  evaluation 
at  the  n level.  This  is  accomplished  by  setting  (3  = 0 in  Eq.  (l4b).  Although 
§ could  be  treated  implicitly  and  linearized,  this  would  unnecessarily  complicate 
the  difference  equations.  The  viscous  terms  involving  V *U  contain  mixed  deriva- 
tives whose  treatment  by  ADI  methods  is  somewhat  awkward,  as  mentioned  previously. 
For  the  solutions  presented  here  and  other  test  cases  computed  while  developing 
the  method,  the  explicit  treatment  of  the  aforementioned  viscous  and  dissipation 
terms  had  no  observable  adverse  affect  on  stability.  All  solutions  presented 
here  were  made  using  three-point  centered  difference  formulas  [i.e.,  ot  = l/2  in 
Eq.  (19)],  and  since  steady  solutions  were  the  primary  objective,  the  backward 
time  difference  form  (|3  = 1)  was  employed  throughout.  In  each  case,  the  refer- 
ence length,  Lr,  was  taken  as  the  duct  width,  and  values  of  0.73  and  1.4  were 
used  for  Pr  and  y,  respectively. 


The  flow  geometry  and  coordinate  system  are  shown  in  Fig.  1.  Boundary  and 
initial  conditions  are  required  to  complete  the  problem  formulation.  It  is 
assumed  that  the  duct  is  fed  from  a large  stagnant  reservoir  and  exhausts  into 
a constant-pressure  reservoir.  As  an  approximation  to  these  conditions,  the 
isentropic  stagnation  temperature,  T , and  stagnation  pressure,  PQ,  are  specified 
as  upstream  boundary  conditions,  and  the  static  pressure,  Ps,  is  specified  as  a 
downstream  condition.  These  quantities  can  be  written  in  nondimensional  form  as 

To=T+  ("^±)m2w2  (26a) 


P0  = 


y M 


I RT 


T \ X/O-X) 


PT 


(26b) 

(26c) 
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The  u and  v velocity  components  are  small  and  were  neglected  in  the  definition 
of  T . Including  the  u and  v components  in  the  definition  of  T0  would  couple 
all  five  governing  equations  at  the  upstream  boundary  during  the  z-direction 
step  of  the  ADI  procedure.  Unless  u and  v were  treated  explicitly,  this  coupling 
would  preclude,  during  this  ADI  step,  the  use  of  the  more  efficient  solution 
technique  previously  described,  in  which  only  three  equations  are  coupled. 
Implicit  boundary  conditions  are  obtained  by  writing  Eq.  (26)  at  the  (n  + l) 
level  and  linearizing  by  the  same  procedure  employed  for  the  governing  equations 
at  interior  points.  The  specified  downstream  static  pressure,  P , was  an  esti- 
mate of  that  required  to  maintain  an  average  nondimensional  velocity  of  1.0  at 
the  duct  entrance.  Additional  boundary  conditions  are  required  at  the  upstream 
and  downstream  boundaries.  The  relation  D|wn+1  = 0 was  used  at  points  adjacent 
to  the  upstream  boundary.  This  boundary  condition  is  equivalent  to  a linear 
extrapolation  of  wn+-'-  on  the  upstream  boundary  from  values  at  the  two  adjacent 
interior  points  (on  a line  in  the  z direction),  and  will  be  referred  to  sub- 
sequently as  implicit  linear  extrapolation. 

As  the  remaining  upstream  boundary  conditions,  the  normal  derivatives  of 
un+^  and  vn+'1'  are  set  equal  to  zero  using  three-point,  second-order,  one-sided 
difference  approximations.  At  the  downstream  boundary,  implicit  linear  extra- 
polation relations  were  used  for  Tn+\  un+-1',  vn+l,  and  w114--1-,  together  with  the 
static  pressure  relation.  No-slip  conditions  were  specified  on  the  walls  of  the 
duct,  and  adiabatic  conditions  were  imposed  by  setting  normal  derivatives  of 
temperature  to  zero  using  three-point,  one-sided  difference  formulas.  In  addi- 
tion, the  wall  density  was  determined  implicitly  using  a three-point,  one-sided 
difference  approximation  of  the  continuity  equation.  Since  the  flow  is  symmetric 
about  the  horizontal  and  vertical  planes  passing  through  the  duct  centerline, 
solutions  were  computed  for  one  quadrant  of  the  duct,  and  symmetry  conditions 
were  imposed  on  these  planes  of  symmetry. 


Stability  Tests 


Computed  Solutions 


Two  sequences  of  solutions  were  computed  for  Re  = 60  with  a 6 x 6 x 6 grid 
and  using  different  time  steps,  to  explore  the  stability  properties  of  the  method 
To  provide  a frame  of  reference,  stability  numbers  N^-^  and  N^g  are  defined  as 
the  ratio  of  the  actual  time  step  At  to  the  maximum  allowable  time  step  as  deter- 
mined by  the  CFL  and  viscous  stability  limits,  respectively,  for  one  dimensional 
uniform  flow  at  the  reference  velocity  and  Mach  number.  These  stability  numbers 
are  given  by 

NCFL=-ff(l+■ii^,  (27a) 
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H-~  Re(Az)2 


(27b) 


Typically,  conditionally  stable  methods  would  require  N^,_T  , N £ 1 for  stability. 

. . V-/  J?  JJ  r" 

Multidimensional  forms  of  Eq.  (27)  are  more  restrictive;  however,  splitting 
techniques  can  be  employed  (cf.  MacCormack  & Paullay,  Ref.  24)  to  recover  the 
one-dimensional  forms.  In  actual  computations  with  a conditionally  stable  method, 
the  stability  limits  would  vary  from  point  to  point  with  the  local  velocity, 

Mach  number,  and  mesh  size;  however,  Eq.  (27)  provides  convenient  reference 
quantities. 


The  initial  conditions  used  for  the  stability  tests  are  those  appropriate  for 
uniform  flow  in  the  z direction  at  the  reference  velocity  and  with  the  specified 
stagnation  pressure  and  temperature.  At  t = 0,  the  no-slip  conditions  were 
applied  and  the  downstream  static  pressure  was  impulsively  reduced  to  P . The 
duct  geometry  has  xp  = y^  = 1,  Zp  = 0.5. 

The  downstream  centerline  value  of  (wn+"*‘  - wn)/At  is  a sensitive  indicator 
of  steady-state  conditions  and  is  shown  in  Fig.  2 for  a sequence  of  solutions 
with  M = 0.44.  It  can  be  seen  that  the  method  gave  stable  solutions  for  test 
cases  in  the  range  N^pp  £ 43.2,  £ 4.4,  and  that  steady  conditions  were  reached 

with  significantly  fewer  time  steps  for  higher  N^pp.  A second  sequence  of  solu- 
tions was  computed  for  M = 0.044,  and  similar  results  were  obtained  (Fig.  3)  for 

hcfl  2 l471<  "t  a 20-6- 

The  transient  accuracy  of  the  M = 0.44  test  cases  can  be  assessed  in  Fig.  4, 
although  as  mentioned  earlier,  transient  accuracy  was  not  an  objective  of  these 
calculations.  The  solutions  in  Fig.  4 are  not  independent  of  the  time  step  (N^pp), 
and  thus  there  is  temporal  truncation  error  in  these  solutions.  However,  curves 
a & b display  far  less  dependence  on  N^pp  than  the  remaining  curves,  and  thus  it 
appears  that  convergence  is  beginning  for  T £2.2.  Steady-state  conditions  are 
reached  in  a nondimensional  time  on  the  order  of  4 for  curves  a & b.  Clearly, 
the  temporal  truncation  error  is  much  greater  for  curves  £ & d (10.8  £ N — 21.6), 
although  these  latter  cases  reached  steady  state  in  fewer  time  steps  (Fig.  2). 
Similar  plots  of  the  transient  solution  for  the  M=0.044  test  cases  are  given  in 
Fig.  5. 


It  should  be  emphasized  that  the  degree  of  transient  accuracy  indicated  in 
Figs.  4 and  5 is  not  a general  indication  of  the  accuracy  achievable  at  a given 
^CFL  or  % the  present  method,  since  these  stability  tests  involve  first- 

order  backward  time  differences,  a very  coarse  mesh,  and  impulsive  starting 
conditions . 
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Solutions  attempted  for  WCFL  = 108  (M=0.44)  and  NqFL  = 14,710  (M  = 0.044) 
were  unstable.  Although  the  precise  cause  of  the  instability  is  unknown,  the 
difficulties  did  appear  to  originate  near  the  boundaries.  Since  solutions  to 
Eq.  (17)  in  general  need  exist  only  for  sufficiently  small  At,  one  possible 
explanation  other  than  the  obvious  one  of  a conventional  instability  is  that 
the  system  of  implicit  difference  equations  is  near  singular  for  these  step 
sizes . 

Although  possibly  by  coincidence,  it  is  noted  that  the  breakdown  occurred 
for  time  steps  which  were  larger  than  the  apparent  time  periods  (t  = 4)  of  the 
physical  transients  as  determined  in  Figs.  4 and  5* 

The  effect  of  mesh  size  on  computed  solutions  was  examined  empirically  for 
the  M=0.44  test  case.  Axial  velocity  profiles  are  compared  in  Fig.  6 for  two 
solutions  having  6x6x6  and  11  x 11  x 11  grids  in  the  computed  quadrant  of 
the  duct;  the  two  solutions  are  in  satisfactory  agreement.  The  error  is  larger 
at  the  duct  entrance,  where  the  inlet  conditions  are  somewhat  severe  for  the 
mesh  spacing  used.  In  Fig.  7,  centerline  velocity  and  pressure  are  shown  for 
three  solutions  having  6,  11,  and  21  mesh  points  in  the  z direction.  This  com- 
parison also  reflects  a reasonable  influence  of  mesh  size  on  the  solutions. 

Solution  With  Moving  Wall 


An  additional  solution  was  computed  for  M=0.44,  Re  = 60  but  with  a large 
secondary  flow  caused  by  moving  two  parallel  duct  walls  in  a direction  perpen- 
dicular to  the  axial  flow  direction;  the  wall  speed  is  about  25  percent  of  the 
average  axial  velocity.  Since  the  flow  is  symmetric  about  a horizontal  plane 
midway  between  the  moving  walls,  only  the  lower  half  of  the  flow  field  was  com- 
puted, and  computer-produced  drawings  of  selected  streamlines  for  this  solution 
are  shown  in  Fig.  8.  No  data  or  other  theoretical  studies  are  available  for  com- 
parison in  this  instance  and  the  computed  solutions  are  presented  largely  as  a 
demonstration  of  stability  in  the  presence  of  large  secondary  flows. 


High  Reynolds  Number  Solutions 
Nonuniform-Grid  Transformation 


The  accuracy  of  solutions  computed  with  a given  number  of  grid  points  can 
often  be  improved  by  using  a nonuniform  grid  spacing  to  ensure  that  grid  points 
are  closely  spaced  in  regions  where  the  solution  varies  rapidly  and  widely 
spaced  elsewhere. 

An  analytical  coordinate  transformation  has  been  devised  by  Roberts  (Ref.  25) 
which  is  an  effective  means  of  introducing  a nonuniform  grid  when  the  steep 
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gradients  occur  near  the  computational  boundaries.  If  N grid  points  are  to  be 
used  in  the  range  0 ^ x ^ 1,  and  if  steep  gradients  are  anticipated  in  a region 
of  thickness  a near  x = 0,  then  Roberts'  transformation  T)  (x)  is  given  by 


77(x)=  N-MN-I)  log  (x+Q+l^/og 


(28) 


where  a2  = l/(l  - a).  The  use  of  equally- spaced  points  in  the  transformed 
coordinate  k provides  resolution  of  both  the  overall  region  0 ^ x ^ 1 and  the 
subregion  0 ^ x ^ a.  The  transformation  (28)  was  employed  in  high  Reynolds 
number  solutions  to  provide  increased  resolution  near  the  duct  entrance  and 
in  boundary  layers  on  the  duct  walls.  Values  of  0.1,  0.1,  and  0.25  were  used  as 
cr  for  the  x,  y,  and  z coordinates,  respectively. 


Artificial  Dissipation 


In  computing  solutions  for  high  Re,  it  was  necessary  to  add  a form  of 
artificial  viscosity  or  dissipation  for  the  axial  flow  direction.  Artificial 
dissipation  in  some  form  is  often  useful  in  practical  calculations  to  stabilize 
the  overall  method  when  boundary  conditions  are  treated  inaccurately,  when 
coarse  mesh  spacing  is  used,  or  in  the  presence  of  discontinuities.  [Von  Mises 
(Ref.  26)  has  shown  that  certain  discontinuities  in  solutions  of  the  Navier- 
Stokes  equations  are  possible  despite  the  presence  of  physical  viscosity  and  heat 
conduction  terms.]  The  need  for  artificial  dissipation  arises  in  certain 
instances  when  centered  spatial  difference  approximations  are  used  for  first 
derivative  terms.  The  use  of  artificial  dissipation  is  thus  a matter  of  spatial 
differencing  technique,  and  is  commonly  employed  in  either  explicit  or  inherent 
form  and  in  both  explicit  and  implicit  difference  schemes.  The  particular  form 
described  here  was  adequate  for  present  purposes  but  is  considered  provisional 
and  is  not  recommended  for  general  use,  since  the  formal  accuracy  is  lowered  to 
first  order  for  the  axial  (z)  coordinate  direction. 

The  dissipation  term  used  here  is  based  on  an  observation  (e.g.,  Roache, 

Ref.  27,  p.  162)  that  for  a linear  model  problem  representing  a one-dimensional 
balance  of  convection  and  diffusion  terms,  solutions  obtained  using  central 
differences  for  the  convection  term  are  well  behaved  provided  the  mesh  Reynolds 
number  Re  = |w|&z  Re  is  £ 2,  but  that  qualitative  inaccuracies  (associated  with 
boundary  conditions)  occur  for  Re^z  > 2.  This  suggests  the  use  of  an  artificial 
viscosity  term  of  the  form  ezD^0>  where 

Re  " Re  ^ ~ ReAz>2  (29) 

ReAz<2 


IwlAz 
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to  ensure  that  the  local  effective  mesh  Reynolds  number  is  no  greater  than  two. 
This  dissipation  term  was  added  to  each  of  the  governing  equations,  with  f taken 
as  p for  the  continuity  equation,  u,  v,  w for  the  respective  x,  y,  and  z momentum 
equations,  and  T for  the  energy  equation.  For  scalar  equations,  the  foregoing 
technique  is  equivalent  to  that  developed  by  Spalding  (Ref.  28)  from  an  argument 
not  involving  artificial  viscosity. 

Computed  Solution 


A solution  is  presented  for  M=0.3,  Re  = 600  and  a duct  geometry  for  which 
x = yx  = 1,  z = 5.  An  11  x 11  x 21  grid  was  used,  along  with  a variable  time 
step  with  N up  to  380  and,  up  to  12. 

Computed  axial  velocity  profiles  at  the  duct  entrance  and  exit  are  shown  in 
Fig.  9 for  this  case.  The  subscript  1 denotes  conditions  at  the  duct  entrance. 

The  axial  variation  of  pressure  ratio  and  Mach  number  is  given  in  Fig.  10.  As 
a check  on  the  solution,  two  additional  pressure  curves  are  shown.  One  curve 
represents  the  pressure  ratio  from  one-dimensional  theory  for  adiabatic, 
frictional,  constant-area  flow  of  a perfect  gas  (Shapiro,  Ref.  29).  Using  this 
theory,  the  pressure  ratio  between  two  points  can  be  evaluated  if  the  Mach 
numbers  are  known;  the  average  Mach  number,  MaVg,  (averaged  over  the  cross 
section)  from  the  Navier-Stokes  solution  was  used  for  this  evaluation.  The  second 
curve  assumes  isentropic  flow  and  constant  stagnation  pressure  along  the  center- 
line  of  the  duct  and  was  evaluated  from  the  isentropic  relations  for  a perfect 
gas  using  the  computed  centerline  Mach  number,  M^.  The  general  agreement  of 
pressure  ratio  distributions  seen  in  Fig.  10  is  an  indication  of  internal 
consistency  in  the  computed  solution.  Since  the  average  inlet  Mach  number  is 
only  0.27,  and  since  the  axial  variation  in  density  is  only  about  8 percent, 
compressibility  effects  should  be  relatively  minor  for  this  solution.  The 
computed  pressure  drop  was  therefore  compared  with  the  experimental  measurements 
of  Beavers,  Sparrow,  and  Magnuson  (Ref.  30)  for  incompressible  flow  and  found 
to  be  in  reasonable  agreement,  considering  the  difference  in  M.  The  results 
of  this  comparison  are  shown  in  Fig.  11. 

In  assessing  the  high  Reynolds  number  solution,  the  question  arises  as  to 
whether  the  artificial  viscosity  destroys  the  accuracy  by  changing  the  effective 
Reynolds  number  of  a viscous  calculation.  Here,  it  should  be  emphasized  that  the 
artificial  viscosity  was  used  only  for  the  axial  coordinate  direction,  where 
viscous  terms  are  generally  unimportant;  second-order  accuracy  was  rigorously 
maintained  for  the  two  transverse  directions,  for  which  viscous  stresses  are  large. 
The  magnitude  of  the  artificial  viscosity  terms  in  the  computed  solutions  was 
examined  a posteriori  and  compared  with  other  terms  in  the  equations.  It  was 
found  that  the  artificial  viscosity  terms  were  no  greater  than  about  2 percent 
of  the  largest  term  in  each  equation,  except  at  grid  points  very  near  the  edges 
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of  the  duct  walls  at  the  entrance.  The  specification  of  constant  stagnation 
pressure  and  temperature  at  the  entrance  along  with  the  no-slip  conditions  on  the 
walls  produces  very  large  gradients  in  this  region,  and  the  artificial  viscosity 
terms  were  as  large  as  15  to  20  percent  there.  It  is  believed,  however,  that  the 
accuracy  was  not  seriously  degraded  by  the  artificial  viscosity  except,  of  course, 
locally  near  the  entrance  edges  of  the  walls.  As  a final  check  on  the  solutions, 
the  mass  flow  rate  was  computed  by  integration  of  w over  the  cross  section  at 
each  axial  location  in  the  duct  and  was  found  to  be  constant  to  within  0.4  per- 
cent . 

The  UNIVA.C  1108  run  time  for  the  solutions  presented  here  was  about  3.5  x 
10“^  minutes  per  grid  point  per  time  step,  which  includes  the  use  of  auxiliary 
storage  devices.  Solution  of  the  implicit  difference  equations  was  performed  in 
double  precision.  Convergence  to  a steady  solution  required  from  20  to  100 
time  steps,  depending  on  how  the  time  step  was  chosen. 


Convergence  Acceleration  Tests 

In  the  solutions  presented  to  this  point,  a constant  or  more  or  less 
monotonically  increasing  time  step  was  used,  and  this  time  step  was  selected 
on  the  basis  of  rough  estimates  of  the  physical  time  scales  present  in  the 
transient  evolution  toward  steady  state.  However,  when  transient  accuracy  is 
of  no  concern  (as  is  the  case  for  computing  steady  solutions),  the  time  step 
can  be  regarded  as  an  iteration  parameter  which  can  be  selected  so  as  to  speed 
up  the  convergence  to  a steady  solution.  In  previous  studies  in  which  ADI 
methods  have  been  used  to  solve  iteratively  scalar  elliptic  equations  such  as 
Poisson's  equation,  it  has  been  established  (Ref.  31)  that  a significant  improve- 
ment in  rate  of  convergence  can  be  achieved  by  the  cyclic  use  of  a sequence  of 
acceleration  parameters  (time  steps)  of  differing  magnitude,  rather  than  a single 
parameter.  This  concept  was  explored  to  a limited  extent  for  application  to  the 
Navier-Stokes  system  of  equations.  A sequence  of  solutions  to  a sample  problem 
was  computed  in  which  single  (constant)  time  steps  of  various  sizes  were  used. 

The  sample  problem  consisted  of  two-dimensional  laminar  flow  in  a straight  channel 
with  a length  to  width  ratio  of  10:1,  with  M = 0.5  and  Re  = 300,  and  using  an 
11  x 11  grid.  By  comparing  the  transient  solutions  and  convergence  rates  for  the 
different  time  steps,  it  was  established  that  the  optimum  single  time  step  for 
this  problem  was  about  0.1  times  the  time  required  for  a particle  to  pass  through 
the  channel  at  the  average  flow  velocity,  and  that  about  25  time  steps  were 
required  to  reach  steady  state  in  this  instance.  A sequence  of  time  steps  was 
then  selected  and  used  cyclically  to  determine  whether  any  increase  in  convergence 
rate  could  be  achieved.  For  the  particular  test  case  and  time  step  sequence 
considered,  there  was  only  a minor  improvement  over  the  convergence  rate  for  the 
optimum  single  time  step.  However,  as  shown  in  Ref.  31  for  Laplace's  equation, 
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the  improvement  in  convergence  rate  brought  about  by  parameter  sequences  is 
greatly  enhanced  when  the  spatial  mesh  is  refined.  Since  the  present  tests  were 
run  with  a very  coarse  mesh  (to  minimize  computing  costs),  the  tests  are  incon- 
clusive with  regard  to  the  convergence  rate  of  fine  mesh  calculations,  which  are 
of  considerable  practical  importance.  Finally,  it  is  noted  that  the  optimum  time 
step  corresponded  to  N of  about  100. 
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APPLICATION  TO  TURBULENT  FLOW  IN  A STRAIGHT  DUCT 


Background 

Since  many  practical  flows  at  high  Reynolds  number  are  turbulent,  it  is 
desirable  to  have  the  capability  to  treat  various  turbulence  models  within  the 
present  numerical  framework.  Consequently,  some  test  calculations  were  made  for 
turbulent  flow  in  the  same  rectangular  duct  geometry  used  in  the  previous  laminar 
calculations,  as  a means  for  exploring  the  applicability  of  the  present  method  to 
the  computation  of  turbulent  flows.  The  straight  duct  geometry  is  convenient  for 
such  an  exploratory  study,  not  only  because  of  its  relative  simplicity,  but  also 
because  there  are  some  previous  experimental  and  numerical  results  (e.g.,  Ahmed  & 
Brundrett,  Ref.  32;  Launder  & Ying,  Ref.  33)  available  for  comparison.  For  the 
present  calculations,  a relatively  simple  turbulence  model,  consisting  of  an  eddy 
viscosity  formulation  and  specified  mixing  length  is  employed.  The  shortcomings  of 
this  approach  (e.g.,  the  difficulty  is  selecting  a suitable  mixing  length  distribution 
in  advance)  are  well  known,  and  a number  of  more  advanced  turbulence  models  are 
currently  being  developed(  cf.  Launder  & Spalding,  Ref.  3^)  wherein  various  turbu^ 
lence  quantities  are  determined  as  the  (numerical)  solution  of  one  or  more 
"turbulence  equations".  A number  of  these  more  advanced  turbulence  models  employ 
a turbulent  viscosity  as  a means  of  relating  the  turbulent  stresses  to  the  mean 
flow,  however,  and  the  present  computational  problem  is  therefore  similar  in 
many  respects  to  that  which  would  occur  with  one  of  the  more  advanced  turbulence 
models.  Thus,  alternative  and  less  restrictive  turbulence  models  can  be  incorpo- 
rated into  the  present  computational  framework  in  a future  study. 


Turbulence  Model 

To  account  for  turbulent  transport  processes,  the  governing  equations  are 
time-averaged  in  the  usual  manner  for  turbulent  flows  (e.g.,  Hinze,  Ref.  35)-  The 
dependent  variables  are  represented  as  the  sum  of  a time-averaged  quantity  denoted 
by  an  overbar  (— ) and  an  instantaneous  fluctuating  quantity  denoted  by  a prime  ( ' ) . 

This  process  of  averaging  produces  turbulent  correlations  which  are  conventionally 
termed  Reynolds  stresses. 

The  time  averaging  is  performed  over  a period  of  time  sufficiently  long  to 
remove  the  random  turbulent  fluctuations,  but  not  so  long  as  to  remove  the  time 
dependence  of  the  mean  flow.  Since  the  present  calculations  are  limited  to  subsonic 
flow,  it  is  reasonable  to  neglect  turbulent  fluctuations  in  density.  In  addition, 
the  assumption  that  the  total  enthalpy  is  a constant  E0  was  made,  and  thus,  pressure 
is  eliminated  from  the  momentum  equations  by  means  of  Eq.  (8).  As  indicated  previously, 
it  is  no  longer  necessary  to  solve  the  energy  equation.  After  time  averaging,  the 
Navier-Stokes  equations  can  be  written  in  a form  identical  to  Eqs.  (3a-b),  except 
that  the  force  due  to  viscous  stress  F now  includes  both  laminar  and  turbulent  stresses. 
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In  discussing  the  turbulence  model,  Cartesian  tensor  notation  is  employed,  and  thus, 
the  subscripts  i and  j,  previously  used  to  denote  a discretized  variable,  will  be 
used  to  denote  the  components  of  velocity  u^_  = u,  v,  w,  and  the  space  coordinates 
Xi  = x,  y,  z.  The  components  of  the  force  due  to  viscous  stress  can  thus  be 
written  as 


dxi  ^ e j j)  — ^ ( /DUjUj) 


'J 


(30) 


where  the  rate  of  strain  tensor  e.y  is  given  by 


eij  = ■^‘(duj/dxj  +duj/dxj) 


(31) 


The'  Reynolds  stress  terms  p u^  ul  are  related  to  the  mean  flow  variables  by  means 
of  the  Boussinesq  eddy  viscosity  concept  (cf.  Hinze,  Ref.  35)  5 from  which 


/»uiu,ia-2/xTeij+  ^/>ukuk8  ij 


(32) 


where  pT  is  the  eddy  or  turbulent  viscosity  and  ey  is  the  mean  flow  rate  of  strain 
tensor.  With  the  further  assumption  that  laminar  viscosity  fluctuations  are 
negligible,  Eq.  (30)  becomes 


Fr2£[WT>5i]- 


ax 


(Muk) 


(33) 


It  remains  to  relate  the  eddy  viscosity  to  the  mean  flow  variables.  This 
is  accomplished  by  means  of  a mixing  length  hypothesis  of  the  form 


fjLr  =/>42(2f  :f)l/2 


(34) 


One  possibility  for  specifying  the  mixing  length  i is  provided  by  the  following 
empirical  one-parameter  family  developed  by  McDonald  and  Camara ta  (Ref.  35)  on  the 
basis  of  experimental  evidence  for  two-dimensional  turbulent  boundary  layers: 


tanh(J!I)..2>s 


(35) 


where  6^  is  the  local  boundary  layer  thickness,  K is  the  von  Karman  constant, 
y is  the  distance  from  the  wall,  and  2>  s is  a sublayer  damping  factor  defined  by 


(36) 
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where  P is  the  normal  probability  function,  and 


1/2  P_ 
H- 


Here  t is  the  local  shear  stress,  y + s 23  and  o+  - 8. 


(37) 


Although  the  mixing  length  profile  (35)  of  Mcbonald  and  Camarata  (Ref. . 35)  was 
developed  for  two-dimensional  boundary  layers,  it  can  be  adapted  for  use  in  a 
rectangular  duct  by  assuming,  for  example,  that  y is  the  distance  to  the  nearest 
wall.  In  a fully  developed  duct  flow  6^  may  be  taken  as  the  duct  half -width, 
while  in  the  developing  region  an  approximate  average  value  of  6-^  could  be  calculated 
from  the  solution  as  it  evolves  in  time.  It  remains  to  specify  the  parameter 
( lco/ 6-u ) . For  two-dimensional  equilibrium  turbulent  boundary  layers,  &w/ has  a value 
near  0.09  (McDonald  & Camarata,  Ref.  35);  for  fully  developed  pipe  and  channel  flows 
•flco/ftb  is  known  to  be  approximately  0.l4  (cf.  Schlichting,  Ref.  37).  The  use  of 
^co/6-5  = 0.14  in  Eq.  (35)  provides  a close  approximation  to  the  mixing  length  distri- 
bution for  a pipe  as  given  by  Schlichting  (Ref.  37).  Within  the  framework  of  the 
proposed  mixing  length  model,  Eq.  (36),  an  empirical  formula  is  required  for  the 
variation  of  Xco/fiq,  from  0.09  in  the  boundary  layer  region  near  the  duct  inlet  to 
0.l4  in  the  fully  developed  region  of  the  duct  flow. 


Sufficient  information  has  been  given  to  apply  the  foregoing  turbulence  model 
to  flow  in  a square  duct.  One  physical  shortcoming  of  the  model,  however,  is  that 
it  is  based  upon  an  isotropic  turbulent  viscosity,  and  it  is  known  (cf.  Launder  & 
Ying,  Ref.  33)  that  models  of  this  type  do  not  predict  the  experimentally  observed 
secondary  flows  generated  by  the  turbulence.  The  problem  of  steady  fully-developed 
turbulent  flow  in  a duct  of  square  cross  section  was  considered  by  Launder  and 
Ying  (Ref.  33),  who  modeled  the  Reynolds  stress  terms  which  induce  the  secondary 
motion.  Brundrett  and  Baines  (Ref.  38)  also  showed  that  transverse  gradients  in  the 
Reynolds  stresses  result  in  the  observed  secondary- flow  velocities  in  the  cross- 
sectional  plane. 

Launder  and  Ying  (Ref.  33)  modeled  the  Reynolds  stress  terms  (u'2  - v'-2)  and 
u'  v'  which  appear  in  the  vorticity  equation.  Here,  the  Launder  - Ying  approach 
is  adapted  for  use  with  the  present  velocity-pressure  formulation  of  the  Havier- 
Stokes  equations  and  this  requires  specification  of  the  stresses  u'2,  v'2  and  u'  v* . 
To  simplify  the  analysis,  it  was  assumed  that  the  approximations  employed  by  Launder 
and  Ying  in  modeling  the  Reynolds  stresses  for  fully-developed  flow  are  also  suffic- 
iently accurate  in  the  developing  region  of  a duct  flow.  The  high  Reynolds  number 
analysis  of  Launder  and  Ying  was  based  on  the  Reynolds  stress  transport  equations, 
and  neglected  the  convection  and  diffusion  terms,  and  the  generation  terms  in  the 
equations  for  u'2  , v’-2  and  u*  v* . In  addition,  since  the  dissipation  of  turbulence 
energy  at  high  Reynolds  numbers  is  isotropic,  the  approximation 

fK  ^L=  -2-8  € 

dxk  dxk  3 IJ 
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was  introduced,  where  e is  the  turbulence  energy  dissipation  rate.  With  the  above 
approximations  the  transport  equations  for  uf2,  v'2  and  u'  v*  reduce  to 
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Following  Launder  and  Ying,  the  approximations  for  the  correlations  between 
pressure  and  velocity-gradient  fluctuations  are  modeled  using  the  form  suggested  by 
Hanjalic  and  Launder  (Ref.  39) > which  was  based  on  earlier  work  by  Rotta  (Ref.  4o). 
There  results 
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and 
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where  k is  the  turbulence  kinetic  energy,  P is  the  rate  of  production  of  turbulence 
kinetic  energy,  cq  and  Cg  are  constants,  and  |3  = (6  Cg-2)/l.l,  The  production  term 
Pk  is  given  by  (Ref.  33), 


PK=  -uw 
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(4i) 


and  for  nearly  equilibrium  turbulence,  one  may  assume  that  e is  equal  to  Pq.„ 
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The  use  of  Eqs.  (39-4l)  results  in  the  following  expressions  for  u '2,  v-2 


and  u’  v*  : 


...  2 _ r _h_  ,.'w *iw_  _ c _k_  p 
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where  the  constants  are 


2(6c2~2) 
3~  ll(C|-2C2) 
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The  remaining  Reynolds  stress  terms  u*  w*  and  v*  w’  are  modeled  using 
Eq.  (32),  i.e., 


uV=-  ;vW="1'Tlr2 


(») 


The  effective  turbulence  kinematic  viscosity  v<j,  = p^/p  for  the  Reynolds  stress 
model  is  obtained  from  the  Prandtl-Kolmogorov  formula: 

vt  = c„k,/2iw  (45) 

and  following  Launder  and  Ying,  e is  approximated  by 

e = cDk3/2/?w  (46) 

where  is  a turbulence  length  scale.  Furthermore,  from  Eq.  (4l), 


^7.2 


[(4“) 


'Ti'ax 


( dw  \2 

ax2 


and  thus,  Eqs.  (42)  reduce  to 


Ciw2 


<4*->2+ 


ax 


C5^W 


[(#>2  + (aw 


L'dx, 


d x2 


>*; 


(47) 


(48a) 
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72  = _ c.  )2+  c.i  2 [(aw)2  + (^)2l 

v ',3x2;  csxw  Lvax,  ax2  -1 

(48b) 

I uV  = 

CiW  (ax^W  C = C3C^/CD>  an<^  C5=  C4CI//CD 

(48c) 

where  c = CgCy/cp  f 

and  c^  = c^Cjy/cj)  . 

i 


The  length  scale  must  be  specified  in  the  above  model. 
(Ref.  33)  used  the  form  suggested  by  Buleev  (Ref.  4l). 


Launder  and  Ying 


dn 

J - 1 f deft 

2 J s 

where  the  notation  is  explained  in  Pig.  12. 


(49) 


An  interesting  observation  can  be  made  about  the  Launder-Ying  turbulence  model 
for  the  case  of  equilibrium  turbulence,  i.e.,  e = P^.  It  is  easily  seen  that 
Eqs.  (45-47)  reduce  to  the  form  of  Eq.  (34),  i.e., 


v 


. (50) 


with  the  implied  mixing  length,  i-m,  given  by 


iirr  (51) 

As  stated  by  Launder  and  Ying  (Ref.  33),  the  Buleev  length  scale  4^.,  Eq.  (49), 
provides  a nonrigorous  but  plausible  method  of  accounting  for  the  influence  of  adjacent 
walls  in  the  duct.  On  further  examination,  the  Buleev  length  scale  is  found  to  have 
the  property  that  the  mixing  length  distribution  from  Eq.  (51)  along  a plane  of 
symmetry  is  nearly  the  same  as  that  for  a fully  developed  pipe  flow  (cf,  Ref.  37), 
and  significantly,  it  yields  -Cm  -*  K y near  the  wall  and  a value  of  0.l4  at  the 
center  of  the  duct. 

The  present  adaptation  of  the  Launder-Ying  model,  Eqs.  (48-49),  provides  the 
correlations  between  transverse  velocity  component  fluctuations  required  in  the 
transverse  momentum  equation  force  terms,  Eq.  (30).  In  the  x-j_  - momentum  equation, 


Fi  = "dxT  (^u'2)' 


wJpM- 


ax- 


(p  u'w) 


axj  ^/xei  j 


(52a) 
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and  in  the  momentum  equation 
_ J5  . ~ I T 

In  these  expressions  the  first  two  terms  are  obtained  from  Eqs.  (48)  and  the 
third  term  is  obtained  from  Eq.  (44),  with  the  turbulent  viscosity  specified 
by  Eq.  (34). 

Special  consideration  must  be  given  to  calculation  of  turbulent  flow  in  the 
vicinity  of  walls  in  view  of  the  large  flow  gradients  which  occur  there.  Since  the 
expense  of  simply  increasing  the  number  of  grid  points  near  a wall  may  be  consider- 
able, an  analytical  wall  function  formulation  has  been  employed  in  the  present  study. 
A set  of  three  universal  velocity  profiles  (cf,  Walz,  Ref.  42)  is  employed  in  the 
wall  region,  corresponding  to  the  laminar  sublayer  (y+  ^4),  a transition  region 
(4  < y+  < 26),  and  the  logarithmic  law-of-the-wall  region  (y+  £ 26): 

u+=  y+  for  y+<  4 (53a) 


a^vl2)"^3(/°v'w')  + 


dx 


(2^e2j) 


(52b) 


u+  = 


In  (1  + y+)  + c2  + [(i-cr  c2  a)y+-  c2]e  ay+ 


for  4 < y*<  26 


(53b) 


1/  = C 1 1 ny+  + c2 


for  f>  26 


(53c) 


where 


and 


II 

c> 

c* 

(54) 

y*=  y(you*Re) 

(55) 

u*  = (rw  /p)'/z 

(56) 
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In  Eq.  (54)  u denotes  the  total  velocity  component  parallel  to  the  wall,  u*  is  the 

j nondimens ional  friction  velocity,  Eq.  (56),  and  tw  is  the  wall  shear  stress.  The 

constants  cp,  Cq?  and  a in  Eqs.  (53)  were  taken  as  2.5,  5.1,  and  0.3,  respectively. 

The  finite  difference  form  of  the  velocity  gradient  at  the  grid  point  adjacent  to 
the  wall  point  is  specified  consistent  with  the  appropriate  -universal  profile  given 
above.  In  the  present  formulation  the  assumption  of  constant  shear  stress  in  the 
immediate  vicinity  of  the  wall  is  utilized.  The  specification  of  the  velocity 
gradient  at  the  grid  point  adjacent  to  the  wall  (where  the  velocity  is  known)  is 
equivalent  to  imposing  a "slip"  velocity  at  the  wall  itself.  Hence,  flow  resolution 
in  the  region  very  near  the  wall  is  sacrificed  to  attain  accuracy  in  the  central 
region  of  the  flow  field.  However,  if  accurate  calculations  in  the  wall  region  are 
required  at  a later  date,  refinements  to  the  wall  function  approach  may  he  implemented 
easily  in  the  present  computational  procedure. 

I 

In  the  incorporation  of  the  turbulence  model  into  the  numerical  solution  procedure, 
it  is,  of  course,  recognized  that  the  turbulent  viscosity  ultimately  depends  on 
the  mean  flow  variables.  It  is  impractical  to  attempt  to  linearize  p,m  rigorously, 
however,  as  the  partial  derivatives  of  with  respect  to  the  mean  flow  variables 
are  not  easily  obtained.  Consequently,  for  the  present  calculations,  pp  was  eval- 
uated at  the  n time  level  after  each  time  step  and  lagged  during  the  implicit  compu- 
tation procedure,  i.e.,  p^  was  treated  numerically  as  a given  function  of  the  space 
coordinates . 
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Computed  Solutions 

The  flow  conditions  present  in  the  measurements  of  Ahmed  and  Brundrett  (Ref* 

32)  of  the  mean  flow  properties  in  the  developing  region  of  a square  duct  were 
adopted  as  a suitable  test  case.  The  experimental  configuration  (Ref.  32)  consisted 
of  a 3-625  in.  square  duct  with  a length  of  18  feet.  Ahmed  and  Brundrett  presented 
isovels  (axial  velocity  contours  in  the  cross-sectional  plane)  at  a Reynolds  number, 
Re^  = prUjJ>]y[/ij,r  based  on  the  modified  length  scale  (where  = 3-625  in-), 

of  approximately  I.065  x 105.  Experimental  static  pressure  distributions  (Fig.  3a 
of  Ref.  32)  in  the  axial  direction  were  given  at  Reynolds  numbers,  Re^,  of 
0.59  x 105,  I.O76  x 105,  1.4l  x 105,  and  I.58  x 10^  for  a duct  with  a gradual  entrance. 
The  static  pressure  distribution  for  Rejyj  = I.O76  x 105  from  Ref.  3^  provides  the 
pressure  drop  required  in  the  calculation  procedure  for  a specified  section  of  the 
duct  length. 

As  a preliminary  step  prior  to  consideration  of  the  developing  flow,  the  fully 
developed  flow  near  the  end  of  the  duct  was  computed  to  assess  the  various  compo- 
nents of  the  foregoing  turbulence  model.  The  eddy  viscosity  turbulence  model,  Eqs. 
(33-3^)?  along  with  the  mixing  length  from  either  Eq.  (35)  or  (51)?  was  first 
employed  without  the  Launder-Ying  modifications  (Eq.  42)  to  obtain  predictions 
without  a turbulence -driven  secondary  flow.  Unfortunately,  at  the  time  of  the 
present  writing,  a satisfactory  solution  could  not  be  obtained  with  the  modified 
Launder-Ying  secondary  flow  model  employed  in  the  present  numerical  procedure,  and 
thus  none  of  the  solutions  presented  contain  a turbulence -driven  secondary  flow. 

The  cause  of  the  difficulty  with  the  secondary  flow  turbulence  model  is  presently 
under  investigation  and  is  possibly  in  the  treatment  of  boundary  conditions  or  the 
assumption  of  a quasisteady  turbulence  model,  which  results  in  an  absence  of  time 
derivatives  in  the  turbulence  model. 

For  the  fully  developed  flow  region,  a duct  section  of  length  5Dg  was  chosen 
and  the  pressure  drop  was  taken  as  10.21  Pa  (1.042  mm  R^O).  The  boundary  conditions 
appropriate  for  fully  developed  flow  are  a specified  constant  static  pressure  at  the 
upstream  and  downstream  boundaries.  The  static  pressure  condition  is  satisfied  by 
linearizing  Eq.  (8)  about  tn  using  the  same  procedure  that  was  employed  for  the 
governing  equations.  In  place  of  the  no- slip  conditions  on  the  duct  walls,  the 
wall  function  analysis  discussed  previously  was  used  for  velocity  components  parallel 
to  the  wall.  Otherwise,  the  boundary  conditions  were  identical  to  those  employed 
in  the  laminar  calculations.  Symmetry  conditions  were  employed  at  the  planes  of 
symmetry,  passing  through  the  duct  centerline,  and  the  solution  was  computed  for 
one  quadrant  of  the  duct. 

Solutions  were  computed  for  a duct  of  length  5I>h  using  a 10  x 10  x 6 grid 
with  the  Roberts  transformation  (Ref.  28)  applied  in  the  transverse  coordinate  direc- 
tions to  improve  accuracy  near  the  walls.  It  should  be  noted  that  reducing  the 
Roberts  parameter  a from  0-5  to  0.1  had  a negligible  effect  on  the  computed  mass 
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flow  rate  and  overall  skin  friction.  Also,  integral  momentum  conservation  checks 
performed  using  the  calculated  wall  shear  stress  and  axial  pressure  gradient  showed 
that  the  calculated  axial  momentum  flux  was  conserved  to  within  2.5  percent. 

Axial  velocity  contours  calculated  with  the  mixing  length  Eq.  (51)  are  shown 
in  Fig.  13  along  with  the  measurements  of  Ahmed  and  Brundrett  (Ref.  32)  which  were 
taken  at  the  location  Z/Dp  = 26.2  at  a Reynolds  number  Rep  = 0.881  x 105.  The 
comparison  clearly  shows  the  distortion  of  the  measured  isovels  due  to  convection 
of  high  momentum  fluid  from  the  duct  center  toward  the  corners,  which  is  presumably 
due  to  the  existence  of  secondary  flows  in  the  duct  (Launder  & Ying,  Ref.  32)  not 
modeled  in  the  calculations.  A comparison  of  calculated  and  measured  values  of  skin 
friction  coefficient  shows  some  interesting  features.  The  Preston  tube  measurements 
of  Ahmed  and  Brundrett  (Ref.  32)  and  the  skin  friction  coefficient  deduced  from  the 
measured  pressure  gradient  in  the  fully  developed  region  are  shown  in  Fig.  l4,  along 
with  the  calculations  of  Launder  and  Ying  (Ref.  33).  Here,  the  skin  friction 
coefficient  is  defined  as  Cf  = Tw/G§prUr  ) > where  tw  is  the  average  wall  shear 
stress,  and  the  Reynolds  number  is  given  by  Re-^  = P^^p/p.  The  present  results 
shown  in  Fig.  l4  were  obtained  using  the  mixing  length  models  described  by  Eqs.  (35) 
and  (51)  at  two  different  Reynolds  numbers  (corresponding  to  two  specified  pressure 
drops).  The  difference  between  the  results  using  mixing  length  Eq.  (35)  and  Eq. 

(51)  is  significant,  and  may  be  attributed  to  the  fact  that  the  Buie ev  mixing  length, 
Eq.  (51),  is  considerably  smaller  than  the  present  adaptation  of  the  McDonald- 
Camarata  length,  Eq.  (35)?  in  the  corner  regions  of  the  duct.  These  results  illus- 
trate the  sensitivity  of  the  predictions  to  the  choice  of  length  scale,  and  perhaps 
reflect  the  shortcomings  of  a turbulence  model  which  requires  a priori  specification 
of  a length  scale  for  a noncircular  duct  flow. 
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TRANSIENT  BEHAVIOR  OF  VELOCITY  TIME  DERIVATIVE  AT  DOWNSTREAM  CENTERLINE 

FOR  DIFFERENT  TIME  STEPS 
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NONDIMENSIONAL  TRANSVERSE  DISTANCE,  y/y., 
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EFFECT  OF  MESH  SIZE  ON  COMPUTED  AXIAL  VELOCITY  PROFILES  AT  x/x-]  = 0.5 
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Terms  in  the  governing  equation  (3)  are  arranged  for  application  of  the 
numerical  method  as  follows : 

HT  = (P,  PU,pv,pw,/>T)  (A-l) 

-d(pu)/dx 

-3(/>u2  + /oT/yM2)/dx  + a2(u/Re)/dx2 
'3(jouv)/<3x  + 32(v/Re)/3x2 
' <3(/ouw)/<3x  + d2 (w/Re)/<3x2 
-d(/>uT)/dx  + (1-  y)pT3u/dx  + d2(rr/RePr)/<3x2 


A-l 


R75-9H363-15 


a2w  = 


ST-  ( 


/ 


- dip w)  /dz 

- dipuw) /dz  + a2(u/Re)/az2 
-<3(/>vw)/3z  + a2(v/Re)/az2 

~d(p\NZ  +PT/YMz)/dz  +a2(w/Re)/dz2 

-d(/oWT)/dz  +(l-X)/?T  dw/az  +d2(>T/RePr)/dz2 


I diV  u) 
’ 3Re  ax 


i a(V'U)  i a(v-u) 

’ 3 Re  ay  ’ 3Re  dz 


yi  X-l)  m2  <Z>  /Re 


A- 2 


(A-4) 
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Stokes  equations  is  presented.  The  method  consists  of  a generalized  implicit 
scheme  which  has  been  linearized  by  Taylor  expansion  about  the  solution  at  the 
known  time  level  to  produce  a set  of  coupled  linear  difference  equations  which 
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are  valid  for  a given  time  step.  To  solve  these  difference  equations,  the 
Douglas-Gunn  procedure  for  generating  alternating-direction  implicit  (ADl) 
schemes  as  perturbations  of  fundamental  implicit  difference  schemes  is 
employed.  The  resulting  sequence  of  narrow  block-banded  systems  can  be 
solved  efficiently  by  standard  block-elimination  methods.  The  method  is  a 
one-step  method,  as  opposed  to  a predictor-corrector  method,  and  requires  no 
iteration  to  compute  the  solution  for  a single  time  step.  The  use  of  both 
Second  and  fourth  order  spatial  differencing  is  discussed.  Test  calculations 
are  presented  for  a three-dimensional  application  to  subsonic  flow  in  a 
straight  duct  with  rectangular  cross  section.  Stability  is  demonstrated  for 
time  steps  which  are  orders  of  magnitude  larger  than  the  maximum  allowable 
time  step  for  conditionally  stable  methods  as  determined  by  the  well  known 
CFL  condition.  The  computational  effort  per  time  step  is  discussed  and  is 
very  approximately  only  twice  that  of  most  explicit  methods.  The  accuracy 
of  computed  solutions  is  examined  by  mesh  refinement  and  comparison  with  other 
analytical  and  experimental  results.  Finally,  some  test  calculations  for 
turbulent  flow  were  made  using  a simple  turbulence  model  consisting,  of  an 
eddy  viscosity  and  specified  mixing  length.  The  results  of  these  calcula- 
tions are  discussed. 


SECURITY  CLASSIFICATION  OF  THIS  PAGEfWfion  Date  Entered) 


